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ABSTRACT 


Fourier transform techniques have been the favored methods in the analysis of 
signals and systems. One major drawback of Fourier methods is the difficulty in 
analyzing transient and/or non-stationary behavior. Recent advances in the field of 
wavelet theory show much promise in alleviating these problems. This thesis considers 
the realizations of the wavelet decomposition and reconstruction algorithms for the 
discrete case. The major discussion will involve both the one and two dimensional 
transforms. We also present a multiple-phase development as a second and possibly a 


preferable method for decomposing signals. 
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I. INTRODUCTION 


The recent advances of wavelet theory have spurned much interest in many 
disciplines: earthquake detection, medical EKG’s, and multiple resolution signal 
processing. We are most interested in the last discipline listed, multiple resolution signal 
processing. Using wavelet theory in this field opens several areas of interest. In the 
one-dimensional case, analysis of short-time duration signals such as radar pulses, can 
possibly lead to better discrimination of similar but distinct radar sources. For the two- 
dimensional case, statistical pattern recognition using wavelets can aid in identifying 
skewed or rotated images. Lastly, for both dimensions, wavelet theory can be used in 
data compression, which has major applications particularly in transmitting data on a 
bandwidth-constrained channel. Traditionally, most means to date use Fourier methods 
to attack the problems in the areas of interest discussed previously. However, for 
transient and/or non-stationary phenomena, Fourier techniques give inadequate results in 
the transform domain, even when we modify them to include a dimension of time. Thus, 
wavelets can serve as an alternative to the conventional Fourier transform techniques for 
analyzing such transient and/or non-stationary behavior. 

This thesis develops one and two-dimensional invertible discrete wavelet transforms 
(DWT). 386MATLAB and PROMATLAB Version 3.5 are the software tools used for 
386 or better PC’s, and for the UNIX-based SUN workstations. We did not use Fortran 


and C programming languages to maximize the portability of the routines, but such codes 


could be developed very easily from this thesis. These routines can be used by other 
researchers for more in-depth analysis of signals and systems. 

Chapter II will cover the basic principles of the wavelet theory, particularly for the 
discrete dyadic case. Although many volumes are dedicated to the field, this cursory 
outline of the theory will provide the basis for notation used throughout the thesis. 
Chapter III will entail some important concepts such as causality and spillover effects that 
must be accounted for in the DWT. We expand these ideas further into the two- | 
dimensional case in Chapter IV. Chapter V proposes another method for decomposing 
the data at the expense of physical memory, the so-called multiple-phase DWT 
(MP/DWT). The MP/DWT possibly approximates the continuous WT, with the property 
of time invariance. Chapter VI gives a narrative description of the MATLAB algorithms 


developed. Finally, we present our conclusions in Chapter VII. 


Il. BASIC WAVELET THEORY 


Since wavelets are a relatively new and less established field to most readers, one 
main question asked is "what are they?" We introduce the topic in here by an example. 

Consider a function f(t) that has discrete levels as shown in Figure 1. We 
decompose f(x) into a lower resolution level. In other words, we desire to smooth (or 
average) out f(x), or to low pass filter the function to a coarser level. Figure 2 depicts 
the resulting decomposition function ,f(x), where « is a scaling factor. For the dyadic 
case, we set a to two. If we compare the original function f(x) to the coarser function 
af(x), we note that the detail is ,d(x) —f(x)-,f(x) as shown in Figure 3. If we had the 
detail and the decomposed signal, we can theoretically reconstruct f(x) without any loss 
of information. The decomposition (and the reconstruction) process can go to any 
arbitrary resolution level m. Usually, we sample the signal at a rate higher than or equal 
to the Nyquist rate, and set the resulting sampled data to be at the highest resolution level 
(say level m=0). All other resolution levels would be in the negative direction. This 
process could be repeated for the next lower resolution level by operating on the current 


level in an iterated manner. 


о 2 a e з 
Figure 1. Example Function f(x) 
Figure 2. Blurred ,f(x) Figure 3. Detail а(х) 


The primary goal is to find a set of orthonormal basis functions that will do the 
smoothing of the original level and retain the details for reconstruction in a more 
practical manner. The low pass filtering process uses a function called the scaling 
function ¢(x) while the detail basis function uses the wavelet basis function (x), both 
of which are also orthogonal to each other. 

The described method is a coarse level in the understanding of the wavelet theory. 
We must now go to a higher resolution level of knowledge to better comprehend the 


following chapters. 


A. THE SCALING FUNCTION ф(х) 

The region L?(R) depicts a vector space where a function f(x) resides and satisfies 
the condition of having finite energy, or the inner product of f(x) with its complex 
conjugate is finite. Consider a family of embedded closed subspaces, Vm, such that the 


infinitely largest subspace is L’(R) while the infinitely smallest subspace {0}. 


+ оо 


ЈР одах = (Јо),РО) ) < = 04) 


(lower resolution ... C V, C V, C V,, C ... higher resolution) 


Additionally, for each vector space V,, in V,,,,, let Wn be an orthogonal complementary 


We show the spaces graphically in 


ү ӨМ... 


Ул such that Vs 
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Figure 4. 





Vector Subspaces 


Figure 4. 


Note that if given V,, and W,, to W,,,,, we can obtain the vector space V4,,,, as shown 
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Figure 5. 


A function f(x) will lie in a vector space V,, if and only if it can be written as: 
f(x) = 6,2? eQ"x-n) (2) 
neZ 

where we adjust the width of ¢(x) for unit grid spacing at resolution level m=0. The 
inner product of ¢ with an n-shifted version of itself must equal zero 
( <ф(х),ф(х-п)> = 0). 

In order to normalize the scaling function for the corresponding resolution level in 
the L^(R) space, multiply the scaling function by a factor 27? while also dilating in the 
spatial domain. The following two equations simplify the notation for the rest of the 


derivations. 


d, = 4,0) » 2! e Q"x) 9) 


Pnn 7 $,,0) 7 2! é Q"x-n) = 6, (x-27n) (4) 


Let P, be the orthogonal projection of a function in L'(R) onto the vector space 
Va. Note that Equation 2 1s a harmonic series representation of the function projected 
onto the V,, space, with c,,, denoting the weighting coefficients for the basis function. 
This is similar to the evaluation of the fourier coefficients for the fourier series expansion 


with its basis functions. 


Can = (FO), Omn? (5) 


Thus Equation 5 can be viewed as a convolution of the function with the mirrored 


scaling function sampled every 2" seconds. 


Cnn i f (x) i $,, lo; (6) 
The recursion property for the scaling function only depends on two adjacent resolution 


levels m and m4 1. Since V, C V,,,, we have: 


ф, = >. 9, 3 Omk ) DET (7) 


We define the filter coefficients in terms of the inner product of Equation 7 as 


m 
қ - |00 %Фх-ба = 2° (ф,, Ф.) (8) 


We are now able to bring the previous equations into a form known in the wavelet 
field as the fundamental dilation equation. This is a two scale difference equation 
relative to the two adjacent levels. 


фо) = У 2 КЮ Ф@х-Ю (9) 


k 


A finite set of h(n) coefficients gives very desirable traits for selecting wavelets 
with compact support. Major research in the field has been in the study and 
determination of the family of scaling functions that are of compact support and 
orthonormal. With the scaling function being orthonormal, we define the following 


properties which are relative easy to prove [Ref. 1: pp. 689-691]: 


[ Фо) ах = 1 10) 


у МЮ =1 (11) 
К 
Y» sre. (12) 
т 2 
2 2. h(k) h(k-2n) = 5(n) (13) 


Going into the spectral domain, the low pass properties of the scaling function are 


obvious from the next four equations. 


bi) = als) |2) where HH = F heP (14) 
k 
P- н |7) (15) 
р=1 A7 
Y: 9-1 - H0-1 (16) 
k 
НО) Р + 1н) Р 1 (17) 


В. THEWAVELET BASIS FUNCTION y(x) 
In a manner analogous to the development of the scaling function, the wavelet basis 
function y(x) determines the Wn vector space. In a similar fashion to Equation 2, we 


define the detail function as follows: 


d(x) =) d,, 27 y (27x -n) (18) 


nez 


The wavelet basis function is of the form of a constant-Q band pass filter function, 
E an octave in bandwidth as the frequency increases an octave. Since W,, 1S 
orthogonal to V», the inner product of each of the spaces’ basis functions should equal 
zero ( < (x) , d(x) > = 0). This relationship is so intertwined that by knowing a 
valid set of ma S, a set of Yma S can be derived for a particular resolution level. 
Consequently, knowing the h(n) coefficients will determine the appropriate coefficients 
for the wavelet basis functions given as g(n). 


We define the wavelet basis in a similar fashion as the dilation property in Equation 


Wx) - 2 Э g(k) ф(2х-®) (19) 


The {y,,(x-n)} is an orthonormal set spanning W,,. The g(n) coefficients share many of 


the properties of the h(n) sequence with the following exceptions: 
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| ұбдас = У 200 = 0 (20) 
k 


wo = eff] ІН (21) 


Mallat [Ref. 1: p. 681] and Daubechies [Ref. 2: p. 944] chose the g(n)'s in terms 


of the h(n)'s as follows: 


gn RENE Ел) (22) 


C. THE DISCRETE WAVELET TRANSFORM 

Note that Equation 6 shows the lower resolution, or the approximation, signal's 
weighting coefficients as a convolution of the time reversal of the h coefficients with the 
signal f(x) sampled every 2" seconds. Using the recursive properties of the scaling 
function, it is easily shown that the approximation coefficients of the lower resolution 
levels also can be determined from approximation coefficients of the higher level. In 


particular: 


11 


Су 
tl 


mn = SFO) Au) 
| >. ( Pin» Pina? P mX 


-Y (6,0, [офи 
= бить RO, Pm 


(23) 
1 
нс 
k 
1 
-22У7 ЖО c,,,,(2n+k) 
k 
1 
-2? S hn-j)c, () 
j 
where h(n) = h(-n). 
The detail weighting coefficients can similarly be determined as 
1 
d,, = 27 E EN-DO (24) 
j 


Recall that the dma s are the detail coefficients of the orthogonal projection of f(x) onto 
the W,, vector space. The entire collection of {dma}, or (c4) U (d, : m > -M} is 


called the discrete wavelet transform. 


1. Decomposition 
Decomposing to any lower level requires that the L^(R) function be initially 
convolved with $,(-x), and sampled at a unity interval to obtain the coefficients, co(n), 


as shown in Equation 2, and graphically depicted next: 
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E « Con 


Figure 6. Determining m=0 Coefficients 


Once the Cp, coefficients are known, we use Equations 23 and 24 in an 
iterated manner to obtain the set of coefficients for any arbitrary level. For the dyadic 
case, this process is a convolution, a two-times decimation of the result, followed by a 


2" scaling factor. 


(em on 
Р i 
ы m 
== 


Figure 7. Decomposition Algorithm 
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2. Reconstruction 
Recall that V,,,— V,€ W,. The c,,4's reside in the V,, space and the d,,,’s 
lie in the W, space. Thus, a projection of a function onto the next higher level should 
equal the sum of the projections onto Vm and Wn. Let Pn and Qn be the projection 


operators onto the vector space. 


P if = PID QI 


P, affi 7 Ye, (D, 6-2 0701) 
2 
P tft = c, (06, (x-2"n) - 


О = у d, (n), (x-2"n) 
Matching the terms in Equation set 25 and doing several transformations of 
variables, we get the following: 
Cin DEN. ic, CQ) (Qn 21) * d, (0)g(n-21)] (26) 
l 


Equation 26 shows that the reconstruction of the approximation coefficients 
of one level higher can be realized by putting zeros between the lower resolution 
coefficients, convolving with the respective H or G filter, adding and then scaling the 


result by 277, We show a block diagram of this basic algorithm in Figure 8. 
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malan 2 ё 


Figure 8. Reconstruction Algorithm 


For the rest of the discussion, we will concentrate primarily on the 
approximation and detail coefficients. It is the properties of these coefficients that are 
of importance in relation to the DWT process. 

The basic procedure would be to decompose the data array, starting from the 
highest resolution level. Next calculate the coefficients for the next lower level, while 
retaining the details, dan, and using the resulting approximation coefficients to calculate 
the next lower level. Repeat this procedure until reaching the lowest desired level where 


now both the approximation and detail coefficients are retained. 
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III. ONE DIMENSIONAL DWT DEVELOPMENT 


A. DECOMPOSITION 


1. Causality 

Recalling that the decomposition algorithm of Figure 7, a finite data vector 
of the sampled signal, Co, can be decomposed to its lower resolution coefficients, c ,, and 
diae This process requires that the data array be convolved with the respective filters. 
These filters are the time-reversed version of the H and G vector array. At the output 
of each convolution, there is a two times decimation. The 2!” scalar multiplication 
normalizes energy the in the L'(R) domain. 

Consider a finite array cy, of length L. Let Co to be a non-zero set only if 
the indices, n, satisfy 0 € n x L-1. Also set the number of h coefficients equal to an 
even number N. By using Equation 23, the convolution process can be depicted as 


shifting the h's by a factor of two each time. 





Shift by two ~———— ———b- Shift by two 


Figure 9. Shift Factor of 2 for N=4 
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In an effort to make the transform a causal system, we desire to make all the 
non-zero data coefficients have indices greater or equal to zero since negative indices 
would equate to negative time. This reconstruction gives two possible sets of indices for 


the h filter coefficients. 


e {h(-N+1),...,h(-2),h(-1),h(0)} 


We allow the negative indices for the filter since we already know the filter 
coefficients (a-priori data). These two sets for the filter coefficients equate to the two 
possible phases that can be used for the decomposition. For notation purposes, we will 
set the largest value of the indices to equate to the type of phase decomposition, phase-0 
or phase-1. 

By using Mallat and Daubechies selection of the g coefficients in Equation 
22, there is a problem that occurs. While the resulting convolution with the h mirror 
filter gives the values for indices that are greater or equal to zero, the output of the detail 
coefficients give values for negative indices. In other words, given we define the h’s 
with phase-O indices, the resulting g coefficient indices are all greater or equal to one, 


as shown in Equation 23. 
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In order to get both the c,,,’s and d,,,’s to fall into the right hand side, we 
must find another relation for g in terms of h, that satisfies the properties in Equations 


18-21. The following equation will serve the purpose: 
(п) = (-1)""3 h(-N+3-n) у. 


2. | Phase Selection 
The decomposition can proceed with the choice of phase to select. The 
selection of the desired decomposition phase and size of the input vector, L, determines 


the total number of coefficients at the lower resolution level. 








phase-0: |{c°_,,}| = — 
(28) 
phase -1: (охшау = Bra 


3. . Zero-Padding of the Data Array 
The basic algorithm uses the dot product between the h (and g) coefficient 
vector with a set of data coefficients, Cm» S, that are obtained by a sliding window on the 
input coefficients, Cm+ı n vector. The sliding window shifts to the right on a generated 
work vector to simplify the programming. This work vector is a zero-padded version 
of the input coefficient vector. Depending on the length of the input data array, the 
selection of the phase, and the size of the filters, the required zero-padding can be 


determined for the beginning and the end of the sequence. 
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If we select a phase-0 decomposition, we will generate the work array by 
initially padding one zero on the left. Additionally, if the size of the data array is even, 
we append another zero to the end of the sequence. For the phase-1 case, we pad the 
work array with one zero at the end of the sequence only if the original array is odd. 
After determining all these conditions, we finally pad the work array with N-2 zeros both 


in the beginning and the end of the intermediate sequence. 


Initially pad with one zero if the 
Phase-O : the number of coefficients in the 
array is EVEN 





Initially pad with one zero if the 
PU | the number of coefficients 1n the 
array 1s ODD. 


Figure 10. Definition of the Work Array 


In esoteric terms, a phase-1 decomposition for each level would be more 
efficient. Practically, this savings in memory is very insignificant, since the number of 
practical resolution levels is approximately [logj( | c4, | )| . Additionally, this reduction 


in the size of the workspace is balanced by a larger workspace during the reconstruction 
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process. Most importantly, desiring only to decompose the data with only the phase-1 
case for each level may miss some properties of a signal, which will be further discussed 


in Chapter V. 


4. | Energy Determination 

One quick method in assessing if the decomposition routine is working 
without doing the reconstruction is to check the energy of the signal in each resolution 
level. Recall that Equations 2 and 18 are harmonic series representations. By using |. 
Poisson's Summation Formula (PSF), the energy in the function in the spatial domain 15 
equal to the total energy stored in the weighting coefficients of the transform domain. 
This property holds true when the basis functions are orthonormal in the L?(R) space. 
Thus, the energy of a higher resolution level's approximation coefficients is equivalent 
to the sum of the energies of the approximation and detail coefficients one resolution 


level down. 


У. + У Ic «a2 (29) 


By normalizing the energy to the input resolution level, m —0, a quick check 
of all the levels will add confidence to the program accuracy. In addition, this energy 
check shows the numeric truncation errors inherent on the processing unit used. The 
display of the energy will provide further insight on possible compression schemes for 


classes of signals, since the energy may be only concentrated at certain resolution levels. 
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Now as an example, consider a one-dimensional transient signal', with a 
phase-1 decomposition for all the levels and using a Daubechies 2-tap filter (N —4) 


[Ref. 2: p. 980]. Тһе following diagrams show the original signal and the energy plot. 


x10—35 
TO 


О SO 100 150 гоо 250 


Figure 11. Transient Signal 


* Professor Michael A. Morgan, Chairman of the Department of Electrical and 
Computer Engineering at the Naval Postgraduate School, provided the transient data 
included in the MATLAB routines called "et90.m ." This signal is the back-scattered 
transient electromagnetic field from a 10cm long thin-wire illuminated by a double- 
Gaussian impulse, computed using a time-domain integral equation. The other transient 
signal provided by Professor Morgan is "et60.m". The numbers 90 and 60 correspond 
to the angle of incidence that the double-Gaussian impulse impinges onto the wire's axis. 
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N 


Norrmolized 


y of the 








“Ca coefficients 


— 4 


Resolution level 


of the ''d' 





Phose: 





coefficients 


11111111 









ло — &з — 6 — 4 — 2 с 


Feesolutiorn level 


Energy in each Level 


The energy check going up from the lowest resolution is depicted next. Note 


that the Daubechies h filter coefficients are only significant to 10^ places, thus they 


cause insignificant numerical errors in the reconstruction. 


Verification of the energy in each resolution level 


Level Energy in Energy in C Difference 
level m-1 level m 
с+а 

-7 0.000024166621151 0.000024166621151 0.000000000000000 
-6 0.000258409042490 0.000258409042489 0.000000000000000 
-5 0.001707402900115 0.001707402900114 0.000000000000001 
-4 0.015878270678385 0.015878270678372 0.000000000000013 
-3 0.197762381708774 0.197762381708634 0.000000000000140 
-2 0.474398688169256 0.474398688168927 0.000000000000329 
-1 0.936969604304463 0.936969604304102 0.000000000000361 

0 1.000000000000806 1.000000000000000 0.000000000000806 
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5. Approximation and Detail Time/Scale Diagrams 

One interesting result of the DWT for the one dimensional case is the display 
of the energy on the time/scale plane. The scale can easily be related to frequency by 
noting that the highest resolution is the sampling frequency at m=O. As the scale 
decreases one level, we halve the center frequency for this scale. 

Figure 7 shows the decimation of the coefficients by a factor of two. When 
we obtain the output of the decomposition process, the coefficient’s are in a packed 
format, meaning that the appropriate time spacing is not preserved. Thus, 2-1 zeros 
must be padded between each coefficient for a particular resolution level m. 

The following time/frequency diagrams show the transient signal, mentioned 


in the previous section, with the proper spacing of the coefficients’ energy in the 


transformed domain. 
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Figure 13. Mesh Plot for the Energy of the Approximation 
Coefficients to Eight Resolution Levels 
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Figure 16. Contour Display for the Energy of the Detail 
Coefficients to Eight Resolution Levels 


В. RECONSTRUCTION 


Figure 8 and Equation 25 show the realization of the recomposition routine by first 
inserting a zero after each coefficient. Next we convolve the array with their respective 
filter, sum, and finally normalize by a factor of 2™?. This process is essentially an 


interpolation. 
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Since we have defined the indexing of the h and g coefficients during the 
decomposition phase, the recomposition process is fixed to one of two methods 
depending on the decomposition phase used to get to that resolution level. The same 
dyadic convolutional algorithm can be used for the recomposition after initially massaging 
the input data arrays, and shifting one unit to the right instead of two. Initially, the h's 
and the g's are reversed in order. 

Now we generate the workspace vectors for the respective coefficients. Append a 
zero AFTER each of the input coefficient vectors. Then add two more zeros at the end 
of the modified work vectors. Finally, if the phase-1 decomposition process was utilized 
to obtain the coefficients, append another zero at the beginning of the modified sequence. 
The diagram below summarizes the work vector definition prior to running the 


convolution algorithm. 


N-2 


Т7 оосо 
pe 420240400! Кох 55% С "n 
25 50506060506 > 555 х 506 OO 
m 422404000! 000000000004 
m po 525252 20095 %% . х e 
RRS XO о 0500000 





Figure 17.  One-dimensional Reconstruction Work Array 
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АП the information needed for the reconstruction is the selected phase 
decomposition to get the c's and d's for that resolution level. One factor to note is that 
with Equations 26, there is a possibility that an extra coefficient be generated in the 
reconstruction process. This extra coefficient will theoretically equal zero, but 
numerically it may not and can be cascaded into a notable error if not accounted for. So 
since we know the length of the input vector at resolution level m —0, the size of the data 
array will be known for all levels since we retain all the detail coefficients. Any extra 
zero-valued coefficient. generated can be identified as such and ignored in going up to 
the next resolution level. 

Errors in the regeneration will be totally manifested as numeric errors of the 
computer processors. The following diagram shows the reconstructed transient signal 


and absolute error with the original. 
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Figure 18. Original Transient Signal and Reconstructed 
versions 
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хто" Absolute Error in the Reconstruction 
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Figure 19. Magnitude of the Absolute Error of the 
Reconstruction. Notice the maximum value. 
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IV. TWO DIMENSIONAL DWT DEVELOPMENT 


A. INTRODUCTION 

Arguments for the DWT can be generated for any higher dimension. Of particular 
concern is the two-dimensional case for image processing. The vector subspaces now 
will reside in the L^(R?) vice L'(R) space. A two-dimensional scaling function $(x, y) 
exists whose scaled dilations and translations for a particular resolution level form an 
orthonormal basis for the vector subspace V,,. By making the two-dimensional scaling 
function separable, (x,y) =¢$(x)¢(y), each vector space can be decomposed as a tensor 


product of two identical subspaces in І/ В). 


и = И И, (30) 


The tensor product can be viewed as the in-phase and quadrature-phase components in 
Ум. Lhe properties depicted in Figures 4 and 5 still apply for the two-dimensional case. 
1 1 
К z Кл ФУ, 


= (Wa D Vn) 8 (Wn v, n 


= (Va È Vn) D (Vn Wn) D (Wn Vn) D (Wn O Wa) 
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Notice that ¢,,, defines a set of orthonormal basis functions in V,,,,! and уши іп №. 
Thus, Equation 31 implies that four sets of equations form the orthonormal basis of V4, 


we summarize in the following table. 


Basis Vector Space Basis Frequency 
Function Coefficient Characteristics 


ФЕ) Ф (У) | Vi © Ул. c, (k,l) Vert Low Pass 
Horiz Low Pass 


Фа Ушу) | Va! G W,! | 08,03) | Vert Low Pass 
Horiz Band Pass 

у(х) а ФУ) | №. O Va! 24 (К) Vert Band Pass 
Horiz Low Pass 

VOO)as Va (Y) | Wa O Wa! 3d, (k,l) Vert Band Pass 
Horiz Band Pass 


Let P,,, !D,, ?D,, and ?D,, denote the projection of a L?(R’) function f(x,y), onto 





vector subspaces V,, © V.,, Vn © Wn,» W, € V,, and W,, € W, respectively. Then in 
a fashion analogous to Equation 25, the vertical and horizontal low pass operation of 
f(x,y) onto resolution level m is as follows: 


РЎ = (FEV) ** фифи) = 2S ^ V c (D) b,(2"x-k) , Q"y-D)) (32) 
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Notice that the weighting coefficients, c,, are now indexed in two-dimensions. 
Similarly, the !d,,, *d,,, and °d,, coefficients can be determined which correspond to the 
low horizontal and high vertical, high horizontal and low vertical, and high horizontal 


and high vertical frequency details. 


30 


Since the three two-dimensional wavelets and the one low pass basis functions are 
separable, the two-dimensional case is very easily realized. The weighting coefficients 
can be determined by first decomposing on the rows, then the columns of the array, or 


vice versa. The general algorithm is depicted in Equations 33-36 and Figure 20. 


С Е гу; > c, (k, D h, (1-2j) h, (k-2i) (33) 
ld (i,j) = 2), > c, (k, D) h, (0-2j) g (k-2i) (34) 
4.4.) - 2), > c (k,1) g, (2j) h, (k-2i) (35) 
i) 2); > c, (k,l)g,(-2j) 2, (k-2i) (36) 


m-1 
Cu D E ы у? ) 2 i а, 
Н, С. |2» 2% а. 

H, H, | 2» /2х e 


Figure 20. 2-D Decomposition Routine 
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В. DECOMPOSITION 

We will use the convention in image processing that the upper left corner is the 
starting and the lower right is the ending data point. This way, we extend the idea of 
causality into two dimensions, where the spillover effects will fall to the right and below 
the input array. 

1. Decomposition Mask 

By looking at the 'd,,, coefficients’ indices between the left and right side of | 

the equation, generate a two-dimensional mask by forming the outer product of the 


vertical and horizontal filter coefficients. 


H,G, = [8,7 - [Aj (37) 
where [h] 2 [h(-N 4-2) ... h(0) h(1)]. Similarly the masks for H,H,, G,H,, and G,G, also 
can be easily determined. 

We require four masks to decompose successfully an image as suggested in 
Figure 20 to generate the coefficients c,,,, !d,,, ^d, ,, and ?d, ,. H,H,, H,G,, G,H,, and 


G,G, are the corresponding filter masks as shown next: 
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h(-2)b(-2998 1: 0h(-2) -Hh(o)n(-2) 105 2) 
n У h У n v h У 
Һ(-1)Һ(-1) по) а(- 1) ве (-1) 

h У h v h v 
ван) ПОШАО) h (0)h(0) h(1)h(0) 
h v h v h v h v 
2) СЕ) (1) ОЕ (Ш) DTI» 
h У h v h У h v 













h(-2)h(-1) 






Figure 21. Mask H,H, 
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h V 






Һ(0)41-2) h(1)g(-2) 
Һ(0)4(-1) x 

ва) h(0)g(0) h(1)g(0) 
h v h v h v 

ви) h(0)g(1) h(1)g(1) 
h v h У h v 


h(-2)g(-2) 


h(-2)g(-1) 


h(-2)g40) 


h(-2)g(1) 
h v 










h(-3)g(C 1) 









Figure 22. Mask H,G, 
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g(-2)H(-2) |[g(-2)h(-2) oe eZ) SIRE 
h v h v h у у 

gí(-2)h(-7) g(-1)h(-1) euo» h(-1) ант т) 
h v У h v h v 


g(-2)h(0) g(-1)h(0) g (0)h(0) g (1)h(0) 
h v h У h У h v 
g(-1)h(1) g (0) h (1) g(1)h(1) 
h v v h v 


Figure 23. Mask G,H, 


9(-2)9(-2) 9(0)9(-2) 9(1)9(-2) 
g(-2)g(-1) g(-1)g(1) g(0)g (-1) g (1)g(-1) 





- - О 
g, ( 2)g(0) gt 1)g(0) ӘЛЕК ) g(1)g(0) 
а(-2)4(0) g(-1)9(1) g(0)g(1) а(1)4(1) 
h v h v h У h v 


Figure 24. Mask G,G, 
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These masks generate the corresponding coefficients by shifting by a factor 
of two to the right and downward through a workspace array. With the one-dimensional 
case, there were only two possible phase decompositions. Now four possible phases 


exist. Obviously, the workspace array is different for each phase decomposition. 





Phase-11 


Figure 25. The Four Possible 2-D Decomposition Phases 
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2. — Lero-padding of the Data Array 

As stated from the previous section, the work space arrays are different from 
each other depending on one of tne four decomposition phases selected, the size of the 
mask, and the size of the input array. These considerations must be noted. 

The generation of the work array for the phase-00 case first consists of 
determining if the input array rows and/or columns are odd or even. If either of them 
are even, append a corresponding extra zeros row or column to the bottom or right of 
the input array, respectively. Given the number of rows and columns for the mask as 
N, and N,, add N,-1 zeros rows above and N,-2 zeros rows below the modified input 
array. Append N,-1 zeros rows to the left and N,-2 zeros rows to the right of the 


modified input array. 



















N -2 
h 
— 
m оіоісіоісісісісіо If # columns is even 
N " 1 MASK] loiolololololololo] add one more column 
TT Mello of zeros 
B b- Jm 
Soo Ср 66 
ооо [11111 оо 
- 2 
V 
N, - ] If # rows is even 
add one more row 
of zeros 


Figure 26. Phase-00 Work Array 
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For the phase-10 case, if the original array's rows are even and/or if the 
columns are odd, add one respective zeros row/column to the corresponding bottom/right 
of the modified array. Next add N,-1 zeros rows above and N,-2 zeros rows below the 


work array. Lastly, add N,-2 zeros columns to the right and left of the array. 












N и - 2 
№ -1 ео 
V of zeros 
АТА | | РО 
ojl | | | tt | Bolo 
||| РОО 
oll в -4--- 
ош ОО 
ооо ооо } N -2 
оо ооо ок о[оіо/оо У 
N, $ 2 If # rows is even 


add one more row 
of zeros 


Figure 27. Phase-10 Work Array 


Alternatively the phase-O1 case, if the original array’s rows are odd and/or 
if the columns are even, add one respective zeros row/column to the corresponding 
bottom/right of the modified array. Next add N,-2 zeros rows above and below the work 
array. Lastly, add N,-1 zeros columns to the right and N,-2 zero columns to the left of 


the array. 
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Figure 28. Phase-01 Work Array 


Lastly for the phase-11 case, if either of the input array's rows or columns 
are odd, then include a corresponding extra zeros row or column to the bottom or right 
of the input array, respectively. Given the number of rows and columns for the mask 
as N, and N,, add N,-2 zeros rows above and below the modified input array. Add N,-2 


zeros rows to the left and to the right of the modified input array. 
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If # columns is odd 
add one more column 


N-21 





of zeros 
-2 
N h 7 2 If 4 rows is odd 
add one more row 
of zeros 


Figure 29. Phase-11 Work Array 


3. | Energy Determination 
Similar to the one-dimensional case, adding up the squares of the four sets 
of coefficients at resolution level m will equal the energy in the coefficients for the next 
higher level, m+1. This is a natural two-dimensional extension of Poisson’s Summation 


Formula, discussed previously in Chapter III.A.5, as shown as follows: 


ЕЕ П Де рар) (38) 
j k j k 
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As an example, consider an image of a centered square, in Figure 30, 
decomposed down to three resolution levels with a phase-00 HAAR wavelet. The energy 


plot in Figure 31 shows the distribution of the c's and the four d's for each level. 


Square Image 3D display of the Image 
100 





50 100 


Figure 30. Square Image 
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Figure 31. Energy Distribution 
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Except for numerical computations of the computer processor, all the energy 
from the lower d coefficients plus the energy of the lowest c coefficients add up to the 
energy of the original coefficients, cy. The following energy check results for the two- 
dimensional case, normalized to cp. Notice that there are no errors since we used Haar 
h filter coefficients and that the two-dimensional decomposition uses a normalization 


factor of two vice the square root of two. 


Verification of the energy in each resolution level 


Level Energy in Energy in C Difference 
m level m-1 level m 
с+а1+а2+а3 
-2 0.440942796610169 0.440942796610169 0.000000000000000 
=1 в. 751059322033898 0.751059322033898 0.000000000000000 
0 1.000000000000000 1.000000000000000 0.000000000000000 


C. RECONSTRUCTION 

Since the scaling and wavelet basis functions are separable, the reconstruction 
algorithm closely follows the same process as in the one-dimensional case. First each 
row/column is zero padded appropriately and reconstructed with the one-dimensional 
technique. Then process each resulting column/row in a similar fashion but with a one 
unit shift to the right and then downward, as in a raster scan. We show this process in 


Figure 32. 
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T у Said 
| 


Figure 32. 2-D Reconstruction 


By taking another look at Figure 32, a more efficient method can be realized. By 
padding each coefficient in the array with one zero in both the horizontal and vertical 
directions, we can then slide a reconstruction mask through the padded array. 

1. Reconstruction Mask 

We generate the reconstruction masks in a similar fashion as in the 
decomposition case, but with the important note that the coefficients are in the reverse 
order to accommodate the reconstruction convolution depicted in Figure 32 and the 


following equation in matrix notation 


GH, = [Á] · [8] (39) 


where h=[h(1) h(0) ... h(-N+2)]. 
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Another simpler way to generate these masks is to reverse the order of the 


decomposition mask in both the horizontal and vertical directions. 


We show the 


comparison of the two masks in Figures 33 and 34 for the G,H, case only, all other 


reconstruction masks can be similarly related. 


90-2) (0-2) ааа) 9 (0)һ{-2) ө (2) һ(-2) 
g(-1)hí(-1) g (0)h(-1) g(1»)n(-1) 
g(-2)h(0) g(-1)h(0) g (0)h(0) с (1) Һ (0) 
һ Vv һ № n v hn v 
© (0)һ (1) qi-1)ncr) g (1)n(1) 
v h м һ v 


Figure 33. Decomposition Mask G,H, 
Suh» 9 (0) Һ(1) g(-1)h(1) 
h У һ v n v 
g(-2)h(0) 
h v 
g(-2)h(-2) 
h м 












g(-2)h(-1) 















с (1) Һ (0) g (0)n(0) g(-1)h(0) 
n v h v n У“ 

g (1)h(-1) g (0)h(-1) g(-1)h(-1) 
h м n v h v 

Sg (1)h(-2) p gf{-1)h(-2) 


Figure 34. Reconstruction Mask H,G, 
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2. — Zero-padding of the Coefficient Arrays 

In order to use the reconstruction masks, we generate the work array for the 
lower level coefficients to accommodate for the decomposed coefficient array and mask 
size, and the decomposition phase that was selected. 

The phase-00 reconstruction work array is the smallest of the four work 
spaces. The coefficient array is padded initially with one zero to the right and bottom 
of each of the individual coefficients. This makes the modified array have an even 
number of rows and columns. Then append N,-2 columns of zeros to the right of the 


array followed by N,-2 rows of zeros on the bottom. 


0 
0 
0 
0 
0 
0 
0 





0 


N -2 columns of zeros 
V 


- Initial Array Coefficients 1n a packed format 





- Padded zero put to the right, bottom or the bottom right of the 
individual packed coefficients 


Figure 35. Phase-00 Reconstruction Work Array 
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For the phase-10 case, add one extra column of zeros to the left of the phase- 
Q0 work array. Alternatively for the phase-01 case, append one extra row of zeros to 
the top of the phase-00 work space. Lastly for the phase-11 situation, insert both one 


row and one column of zeros to the top and left of the phase-00 case. 


Add one row of zeros for 
either: 


Phase-O1 
Phase-11 


Phase-OO 
Reconstruction 
Work Array 


from the previous 
diagram 





Add one column of zeros for 
either: 

Phase-10 

Phase-11 


Figure 36. Work Space for Phase-10,10,11 Cases 
3. | Two Dimensional Example 
Consider the diagram of Figure 30, the following plots are the resulting mesh 
display of the decomposed coefficients three resolution levels down (m —-3) using a Haar 
wavelet and a constant phase-00 selection for each level. The top left diagram is the 
approximation coefficients, Cm. Going in a clockwise manner, the following plots equate 


to !d,, ^d, and ?*d,, coefficients. 
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Figure 37. Mesh Display of the Coefficients for level m=-3 
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Figure 38. Contour Display of the Coefficients for level 
m--3 
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Figures 39 and 40 compare the original and reconstructed image coefficients 


displayed as mesh and contour plots. 
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Figure 39. Reconstructed and Original Mesh Comparisons 
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Figure 40. Reconstructed and Original Contour Comparisons 
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V. MULTIPLE PHASE DWT 


For both the one and two-dimensional transforms, we selected a single phase for 
decomposing to the next lower resolution level. If we decompose down to a fixed level 
m-M while picking different decomposition phases for each level, we generate a 
different set of coefficients. This set of decomposition coefficients and any other set of 
decomposition coefficients will successfully reconstruct the original input array as long 
as we note the phase used to obtain each particular resolution level from the previous 
one. Such a disparity in the sets of coefficients suggests that the DWT process of a two 
times decimation is not shift-invariant. In other words, if we decomposed the input array 
with one type of decomposition scheme, our coefficients would be different if the input 
array is shifted before the low pass and band pass filtering processes occurs. As a note, 
using the energy check routine, we can see the distribution of the energy of the 
coefficients change as the phase decomposition scheme varies. So this shift variant 
property can be further exploited by varying the phase for possibly optimizing or 
concentrating the energy of families of signals to certain resolution levels. 

However, we desire linear shift invariant (LSI) systems for the analysis of data. 
For the continuous WT, this property is generally true, however for the DWT case, this 
aspect is not. One way to bypass this obstacle is to account for all the phases during 


each decomposition. We call this technique the multiple-phase DWT (MP/DWT). A 
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thorough discussion for the one-dimensional case will first be presented, which will then 


be extended into the two-dimensional case. 


A. ONE-DIMENSIONAL MP/DWT 

Consider the case of going from the data input resolution m=0 to m=-1. The 
coefficients that can be generated are one of two sets, depending on whether we selected 
the рһазе-0 ог the phase-1 case. We really need only one of the two sets to decompose 
to the next lower resolution level m=-2 to maintain perfect reconstruction. Also, we 
need to pick which phase to decompose from m=-1 to m=-2. Thus for level m=-2, 
there are four possible phase combinations that can result: 00, 10, 01, and 11. The 
ordering of the subsequent phases is read from left to right to correspond to decomposing 
from level m=0 to m=-2 and we call this a phase vector. Thus пос. would be 
interpreted as the set of approximation coefficients at resolution level minus two and 
decomposed with a phase vector of [10] with an index of five. In general, for any 
resolution level m <Q, the total number of phase combinations is 2™. 

Looking again at going from the m=0 to m=-1 case only, we can see that by 
sliding the filter coefficients by one instead of two units to the right, the resulting 
coefficients alternate on the phase selection starting with the first coefficient in the phase- 
0 set followed by the first coefficient in the phase-1 set as shown on the next diagram. 
We determine the total number of coefficients in level m=-1 just like a linear discrete 
convolution, which is the number of data points from level m=0 plus the length of the 


filter mask minus one (| {co} | + N- 1). 
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Figure 41. Alternating of Coefficients by Phase 


Going from level m=-1 to m=-2 will require a slight modification before shifting 
the filter mask by one to the right as discussed in the previous paragraph. Since two 
sets of decomposition coefficients are interleaved, a zero must be appended between each 
coefficient in the filter mask for properly decomposing each set of coefficients. The first 
four values in the resulting set correspond to the phases: 00, 10, 01 and 11. This 
association repeats for each four values in level m=-2. In fact, the spacing of each 
coefficient in a desired phase is properly indexed by the same amount as the number of 
zeros padded for display purposes in the Chapter III. A.5. 

This observation can be generalized to any arbitrary level m. In this case, separate 
the filter coefficients by 2?-1 zeros before the conducting single shift to the right. By 
noting the resolution level m, the order of the phases can be determined as the binary bit 
reversal of the binary values counting from zero to 27-1. We show this characteristic 


as a phase-tree diagram for resolution levels m=0 to m=-3. 
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Resolution Level 


© mc 
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1 пі --1 


m = -2 





Figure 42. Phase-Tree Diagram for Resolution Levels m=0 to 
m--3 


We now define new filter sequences, h,, and g,, (m <0) for decomposing the values 


Cm, to the next level c,,,,. The revised filter sequences are defined as 


h, = [ h(-N+2) (0), h(-N+1) ©}... ©} A(O) (0), A(1) ] т 


2, = [ £(-N*2) (0), g(-N*1) (0), ... (0), g(0) (0), e(1) ] 


where {0},,=2™-1 length zeros vector 
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Thus, the MP/DWT equations are given as follows [Ref. 3: p. 3][Ref. 4]: 


1 


Cn-li B |» V2h, (K) Сы (41) 
к--М»2-(М-1)27"-1) 
and ; 
d, -l,i = М? 28 KK д -k-1 (42) 


k=-N+2-(N-1)(2™-1) 


The total number of coefficients also can be generalized in terms of the size of the 
input array, | dy | , the size of the h or g mirror filters, N, and the current resolution 


level, m: 
lent = leg] + W-1)2™-1) (43) 


So enough memory must be allocated if we desire lower and lower resolutions, since the 
number of possible phases increases by a power of two each step downward. 

Consider a BPSK signal shown in Figure 43. Using a Haar wavelet and a phase 
vector of [11111111], a contour display of the energy in the approximation coefficients 
and of the detail coefficients results and are shown next. Notice that much of the 180° 


phase shifts are not readily detectable. 
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BPSK Signal 


O =o 100 1900 200 zc 
Figure 43. Sample BPSK Signal (Ref. 4] 
zoomed Contour Display of the Approximation Prose: 11111114 
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Figure 44. Approximation |c 
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Zoomed Contour Display of the Detoil Prose: 11111111 
з 


= resolution level 


0 1 QQ 150 zoo 250 
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Figure 45. Detail |d,,|’ 


Now for the MP/DWT case, the 180° phase shifts of the BPSK signal are now very 


apparent in both the approximation and detail coefficients: 


Multiphnose Contour Displey of the Approximotion 


-n resolution level 





10 contour levels 


Figure 46. MP/DWT |с, |? 54 


Multiphose Contour Display of the Detail 


-n resolution level 
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то contour levels 
Figure 47.  MP/DWT |d,|? 

Twice the energy of the coefficients for the current resolution level will equal the 
total energy on adjacent lower level's approximation and detail coefficients, since there 
are two valid sets of coefficients for the reconstruction to the higher resolution. So the 
energy checking routine discussed in the one-dimensional DWT is still valid as long as 
we account for a factor of one half when going up one resolution level. 

The actual zero padding of the filter masks is not the most optimum method in the 
determination of the decomposition coefficients’ proper phase order. A much more 
efficient technique can be determined by noting the pattern that develops as Equations 
42 and 43 and writing it out explicitly. The following diagram shows the developing 


pattern for m=-2. 
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Са Са 

с(0) -а3с(0) 0 0 

fe =a3c(1) 0 0 

Gta) =a3c(2)+a2c(0) 0 

= 05) =а3с(3)-+а2Е(1) 0 

ССЗ) =a3c(4)+a2c(2)+aic(0) 

CS.) =а зе (ы у айе (БББ айс г 

cio) =a3c(6)+a2c(4)+a1c(2)+a0c(0) where: 


C) =а3с(7) +а2с(5) +а1с(3) +абс(1 ) 
[a0 a1 a2 a3 ... aN-1] 


C(L*N-3) =a3c(L-1)+a2c(L-3)+a1c(L-5)+a0c(L-7) [h(-N+2) mento) rei 
or 


C(L*N-2) 0 a2c(L-2)*a4c(L-4)*«a0c(L-6) 
C(L*N-1) 0 а2с(1-1) +а1с(1-3) +абс(1.-5) 
c(L+N) 0 0 dadic(b-2)*a0c05-4) 
c(L+2N-3) 0 0 alccb-1) 0905853) 
EIE SZ NS 0 0 0 абе С Е=4 ) 
СТАНЕ 0 0 0 ӨШЕСсіІ:-1) 


жк 


[8(-№+2) ... (0) О 





Figure 48. MP/DWT Pattern Development for m=-2 


Notice that the coefficient set, c,,,, is a summation of N vectors with | c | --(N- 
1)(2™-1) elements. We pad some of the vectors with zeros either in the beginning or the 
ending of the approximation sequence of the current level and multiply by an appropriate 
scalar value from one of the filter coefficients. The following equations summarizes the 


development: 


{Cm-1} = а, Г, а-ы ү2%(-Мз2зр) (44) 
p=0 
N-1 

а i] т bl w= y29(-N+2+p) (45) 
p=0 
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where 


ү, |04) (6. 0.) 49 


(0,,] — (zeros coefficient row vector of size (N-1-p)*2™} 
(0,.) — (zeros coefficient row vector of size p*2™}. 


This vectorized technique is faster than the iterative process of padding the filter 
coefficients appropriately. However, we need more buffer space to calculate the data, 
which is dependent on the number of resolution levels desired. 

The reconstruction process works in a very similar fashion as in the DWT case. 
All we need is the selection of the coefficient sets, which in turn determines the unique 
path back upwards in the phase-tree diagram, similar to Figure 8. Once we select the 
desired phase at the lowest resolution level, the we extract the corresponding 
approximation and detail coefficients and put them in a packed format. Then the DWT 
reconstruction commences to get up one level, where we extract the next set of detail 


coefficients. We repeat this process until we reach resolution level m=0. 


В. TWO-DIMENSIONAL MP/DWT 

Since we define the basis functions as separable in the two orthogonal directions, 
ме сап и5е (һе same arguments discussed for the one-dimensional case. When we talk 
about coefficients and masks, they are now two-dimensional arrays as discussed in 


Chapter IV. Thus, the basic equations for the two-dimensional MP/DWT are 
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1 
с 20-2 ЖТ H,H ((j,é (a+j-1,b+k-1) 
J7-N,*2-(N,-1)27"-1) 
k- -N,*2-(N,-1)27"-1) 


(47) 


1 
‹ n CDE? у Уу H,G Qk) (a+j-1,b+k-1) 
j*-N,*2-(N,-1)2"-1) 
--М,2-(М,-1)27"-1) 


(48) 


1 
24 (а) =2 УУ С.Н. (,Ю& (а+)-1Ь+К-1) 
j*-N,*2-(N,-1)27" -1) 
k--N,*2-(N,-1)27"-1) 


(49) 


1 
"ge (ab) уу G,G .(,)€, (a *j -1,b *k- 1) 
Jj*-N,*2-(N,-1)2 ^ -1) 
--М,я2-(М,-1)27"-1) 


(50) 


where the top left corner of the mask corresponds to location (a,b) —(0,0), and the 
positive directions are to the right and downward [Ref. 3: pp. 2-3]:. 

Even the one-dimensional vectorized process converts into a "matricized" process. 
Here, the resulting coefficient array is the matrix addition of N,xN, total matrices, each 


appropriately padded with zeros on all four sides and multiplied by one of N,xN, scalar 


values. The two-dimensional equations result as follows: 


(саа = qué a,, = 2H,H,(-N+2+p,-N+2+q) (51) 
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| т-1 ^ у; Б 4” be = 2H,G (-N+2+p,-N+2+q) 
р=0 4-0 
N,-1 N,-1 

E zi Cog 2, Cro ig 2G, H (SN *2*p,-N*2*q) 
p=0 4-0 
N,-1 N,-1 

: m-1 B Pq Z "m É 2G,G (-N+2+p,-N+2+q) 
p=0 4-0 


= 
(М.-1-р)2 Zeroes 


-m 
Р2 zeroes 
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ЖЕТСЕЗ ЛЕВЫЕ 





Figure 49. Definition of Zn 
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-- zero coefficient array 


(52) 


(55) 


(54) 


The phase-tree diagram can be viewed in three dimensions to more appropriately 
show the proper interleaving of the phases. Since there are four possible phases, there 
will be four more locations in the indexing of the data on the next lower resolution level. 


The total number of coefficients would be 


[Rows * (N,- 1)(2" -1] [Cols, *(N, - 1)(2" - 1)] (55) 


where Rows, and Cols, are the number of row and columns of the original 1mage. 
In order to view the phase-tree diagram correctly, we must consider each resolution level 


as a new plane for the ordering of the coefficients. 
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Figure 50. Phase-Tree Diagram 
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The MP energy display of the coefficients closely resembles the display for the 
single phase case. This energy display differs by the fact that they are the normalized 
average energy values for each level. As before, we normalize the energy to the 
resolution level m=O. Consider as an example, a MP/DWT of the image in Chapter IV 
using Haar wavelets and decomposing down to m=-3. The following display shows the 


MP/DWT energy plot: 


wavelet horiz: Hoar wavelet vert: Н оог 
1 а 
о. = о. = 
ге е 
4a . + 
©». > сә. 
© © 
о © 
= c 1 
dz с 5 
1 1 
eS B 
= е 
4 Q.-4 
© oO. 
© 
O © 


тз шір Њо >= cose 


Figure 51. MP/DWT Energy Display 


61 


and the energy check is as follows: 


Verification of the energy in each resolution level 


Level Energy in Energy in C Difference 
m level m-1- level m 
с+а1+а2+а3 
=? 0.441935911016949 0.441935911016949 0.000000000000000 
zd 0.751059322033898 0. 251059322033898 0.000000000000000 
O 1.000000000000000 1.000000000000000 O. 000000000000000 


The MP/DWT of the three-dimensional energy display for all the coefficients of 


the image and the corresponding contour images for the m=-3 case are as follows: 


wavelet horiz: Haar wavelet vert: Haar 





multiphase case Resolution level -3 


Figure 52. Mesh Display of the Stored Energy for m=-3 
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wavelet horiz: Haar wavelet vert: Haar 





50 
multiphase case Resolution level -3 
Figure 53. Contour Display of the Stored Energy for m=-3 
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VI. MATLAB DWT ROUTINE DESCRIPTION 


The Appendix lists the developed DWT algorithms. We categorize them primarily 
into two classifications, the one and two-dimensional cases. In each situation, we further 
define two subgroups, the single phase DWT and the MP/DWT. Included with the 
software is an ASCII version of this chapter in a README.TXT file. This chapter first 
describes the system configuration to run the routines. Then we list brief steps for each 
of the two general cases. These steps are not totally inclusive, but the program's 
prompts wil] fill any other gaps in the routine. Finally, we discuss important 
considerations for obtaining graphical output of the results since this is highly dependent 


on the system hardware. 


А. SYSTEM CONFIGURATION 

Minimum hardware requirements for the PC version are an Intel-386 
microprocessor with a math co-processor, Microsoft DOS or Digital Research DOS 5.0, 
VGA, at least three megabytes of hard disk space, and four megabytes of RAM, for most 
of the DWT applications. However, we recommend at least eight megabytes of RAM 
for using arrays greater than 10,000 elements, or if m<-5 in the two-dimensional 
MP/DWT case. Also, some DOS memory managers may conflict with the Pharlap DOS 
extenders that MATLAB Version 3.5 uses. For the Sun workstations, in order to operate 


PROMATLAB Version 3.5, we recommend the Open Windows or Sunview graphical 


environments with at least eight megabytes of RAM. Any further memory can be 
negotiated with the system administrator. Note that UNIX-based operating systems 


discriminate characters between upper and lower case letters. 


В. GENERAL PROCEDURAL STEPS 

Initially, we must be in the MATLAB environment prior to running any of the 
routines. Either generate in MATLAB or load the sampled data for resolution level 
m=Q, and assign a variable name. If and files with the prefix, "leg", and the suffix, 
" met" (all in lower-case letters), exists in the current directory, they will be deleted once 
the DWT routines commence. Thus if you desire to retain these files, rename them. 


The DWT routines use these files as the output graphic files for the routines. 


1. One-Dimensional General Procedures 

a. If you have no data, the routines include three files as data that can be loaded. 
Type "load" followed by a space and then one of the three names: "etó0.m", 
"et90.m", or "bpsk". The variable name will be either "et60", "et90", "bpsk". 

b. Invoke "wavld" in lower-case letters. 

c. Enter the variable name of the input data. 

d. Enter the sampling frequency, « RET» if not known. 

e. Select the desired h coefficients. 

f. Select "Y" if you desire the MP/DWT routine and skip to step k. 

g. In you did not select the MP/DWT routine, the single phase DWT commences, 
decomposing the coefficients until resolution level m—- [log;( | c, | )] with an 


initial query of if you desire to see the intermediate resolution level displays. 


h. Choose either the phase-0 or phase-1 decomposition for each resolution level. 
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After the final decomposition, the energy checks commence. 


Select "Y" if you desire to do the reconstruction process. If so, answer yes or no 
if you want to display the intermediate levels. 


Lastly, the time/scale display commences, and the series of questions ask if you 
want to zoom in or out and what number of contour levels you desire. 


For the MP/DWT case, steps i to k are essentially the same except for the 


reconstruction, where you must initially select the desired phase, thereby fixing 
the phase path back up to level m=O. 


oe Two-Dimensional Procedures 


If you do not have any input data, invoke "idata" for a small selection of images. 
The variable "1m" (in lower-case letters) will be the input for step c. 


Invoke "wav2d" in lower-case letters. 

Enter the variable name of the input data array. 

Select the desired h coefficients in the horizontal direction. 
Select the desired h coefficients in the vertical direction. 
Select "Y" if you desire the MP/DWT routine. 


Enter the number of resolution levels to go down. An initial good selection is 
four since this will cover all four possible phases. 


For the single phase DWT case, select one of four of the decomposition phases 
for each level. 


After the final decomposition, the energy checks commence. 


For the single phase DWT case, select "Y" if you desire the reconstruction 
routine. 
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C. OBTAINING GRAPHICAL OUTPUT 

Except for the initial time/frequency plots, all other displays will query you if you 
want to store the plot. If so, the plot is stored as a metafile with the prefix of "leg" and 
the suffix ". met". A number starting from zero is put between the prefix and suffix, for 
example, leg14.met would correspond to the fifteenth stored plot in the DWT routine. 
Output of these metafiles are hardware and system dependent. Refer to the GPP 
command in the MATLAB User’s Guide for more detailed information. Two 
C-shell script files, called "metal3" or "metal4", will generate a temporary Postscript 
file for plotting on a Postscript printer all the metafiles retained from the DWT routine. 
The use of these files is currently only guaranteed to work at the US Naval Postgraduate 
School Electrical Engineering Department's Sun workstations, and we list them here as 
a guide in generating other C-shell script files for a particular system. For the DOS 
version, the following command will generate all the metafiles onto a HP Laserjet III 
printer as long as the we invoke the following DOS command in the same directory as 
the metafiles: for %f in (leg*.met) do gpp386 %f /djet150 /ol /fprn . If you put it in 
a DOS batch file, echo the percent sign ( %% instead of % ). Please note that the 
computer may take some time generating these plots. Be also aware that depending on 


the size of the plot, the graphic files may exceed the printer buffer memory. 
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VII. CONCLUSIONS 


In this thesis, we developed discrete wavelet transform algorithms for both the one 
and two-dimensional cases. The iterative mechanism of decomposing the data 
coefficients by a convolution-and-subsample process is achievable. These realizations 
require that the scaling and wavelet functions form orthonormal sets with compact 
support. Also for the two-dimensional case, the basis functions were separable in the 
two orthogonal directions. 

We found that the output of the decomposition coefficients can be causal provided 
that the definitions of the indices for both filter coefficients, h and g, must be selected 
to accommodate the causal condition. In fact, the definition of the g filter coefficients 
must be adjusted from the definitions given by Mallat and Daubechies. We found for 
the one-dimensional case, the decomposition is simply a dot product of the filter 
coefficients with a set of values windowed from the higher resolution level approximation 
coefficients. This window subsequently shifts two units to the right for the next 
decomposed coefficient. We found that the reconstruction process uses the same 
decomposition algorithm if we separate each of the individual input values by a zero, the 
filter coefficients are time-reversed, and the shift is now one unit to the right at a time. 

The two-dimensional case was a natural extension of the one-dimensional process. 
We found that the decomposition can be realized as a mask element multiplication and 


summation, with the mask shifting two units to the right until the row completion and 
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then two units down to process the next row of coefficients. The shifting of the mask 
is similar to a raster scan. The reconstruction follows the same process as in the 
decomposition but each of the input array coefficients are separated by a row and a 
column of zeros, the filter masks are spatially reversed in both direction, and the raster 
shift is one unit instead of two. 

We found that there were two possible phases in decomposing the data in the one- 
dimensional case, and four possible phases for the two-dimensional case. This led to the 
conclusion the current decomposition algorithm was very shift variant, contrary to the 
continuous version to the wavelet transform, which is shift invariant. We then 
determined that the discrete case can possibly approximate the continuous version if we 
account for all the phases during the decomposition. We then developed the multiple- 
phase discrete wavelet transform for both one and two dimensions. The coefficient sets 
for a particular phase can be interleaved with other sets of a different phase. We found 
that the multiple-phase was more capable in detecting discontinuous properties of data 


than in the regular discrete wavelet transform, such as a binary phase shifting function. 
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APPENDIX DWT MATLAB FILE LISTINGS 


A. ONE-DIMENSIONAL ROUTINES 


(6 ФФФФФФФФЖФФФеФФФФФоФеб06еееФтеететеететекееоееквзкоФФЖЕНИ ЕК ФФ же ооо 


% 15 ЗЕР 92 
% WAVID.M _ The One Dimensional Discrete Wavelet Transform 
% 

% 

* By: LT J. E. Legaspi 

% Professor Lam 

% US Naval Postgraduate School, Monterey CA 

% e-mail: legaspi@ece.nps.navy.mil 

% lam@ecce.nps.navy.mil 

% Phone: (408) 646-3044/2772 

% 


ооо р 


% 

* This routine docs a one-dimensional dicrete wavelet decomposition 
% of a one-dimensionaldata array. The algorithm allows for the 

% selection of the desired phase for the respective decompositon 

* level. The following routines are necessary for the proper 

% operation of this algorithm: 


% 

%  chkputer.m -- algorithm checks for compatible hardware 

%  daubdata.m -- determines the h coefficients for a daubechies 
% compactly supported orthogonal wavelet 

* mp.m -- MultiPhase decomposition 

Ф phsl.m  -- selects the phase-1 decomposition 

* phsO.m  --selects the phase-0 decomposition 

*  dpscoef.m -- display the coefficients for each level 

% enrgld.m -- determines the energy check for each resolution 
Ф reconld.m -- reconstruction algorithm for verification 

% plane.m  -- displays the coefficients in the time-scale domain 
% strpli.m -- stores the plots in meta-files 
j—————————————— 

clc 


chkputer % Checking for the basic hardware necessary 


ое LS ee ОНЕР ОИНИ 

% с0 corresponds to the resolution level 0, the original data 
% array. 

cle 


c0 z input(' Enter the name of the data in vector row format: °); 


% We will check if the data is in row format, if in column format, 
% take the transpose of the input. 


[il jl] =size(cO); 
if jl==1, c0=c0’; 
elseif 11>14)1>1, 
disp('Bad Input Data, Restart the Program’) 
return 
end 
clear il jl 
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fsamp — input('Enter the sampling frequency (Hz): °); 


* The following code determines the desired wavelet 


for 1=1:8,disp((’ ’}),end 
ql 7 ['Enter the number for the desired choice: ° 


*(1) Haar h coefficients : 

'(2) Daubechies h coefficients : 

'(3) Own set of h coefficients : 

' 
disp(q1) 


qqq=input(’ ’); 


if qqq==3, 
h=input(’ Enter your own set of h coefficients’); 
% Checking if the data are compact support 
a=sum(h);b=h*h’; 
1011 = 1е-12; 
while abs(1-8) » toll & abs(.5-b) » toll 
disp(" Your h coefficients are NOT compactly supported!') 
h=input(’ Re-enter the h coefficients or ctrl-z to stop’); 
end 
wwsvlet Z'OTHER'; 
elseif qqq==2, 
qq=input(’How many taps (2-10)? ’); 
h =daubdata(qq); 
wwavelet- ['Daubechies ' ,num2str(qq), '-tap']; 
ele hz2[.5.5];wwavelet Z'HAAR'; 
end 


h=sqrt(2)*h; 9 The factor of sqrt(2) normalizes the energy 


ть 
% Some definitions of variables used throughout the entire 
% algorithm: 

% 


% wwavelet -- name of the selected wavelet 

* LL -- the length of the input array 

% Nh -- the number of h and thus g coefficients 

* lowest -- the number of resolution levels below the input 

% pltcnt -- plot counter for storing desired plots 

%Һ -- the "h" coefficients 

%% -- the *g* coefficients 

% lowest -- the lowest resolution level 

% phsvct -- the phase vector for recording the appropriate phase 
% zro -- Zero vector to record the appropriate zero padding 
% for each resolution level. 


% We will intialize some of the constants: 


LL=length(c0); 
Nh=length(h); 

lowest = -ceil(log(LL)/log(2)); 
№М=0; 

pltent 20; 


Т1 


* Generating the g coefficients 


g 7 fliplr(h); 

for i=2:2:Nh 
g(1.i) = -g(1,i); 

end 


% The following set of lines branches to the multiphase 
% decomposition algorithm if desired. Program ends after 
% the multiphase execution 
cle 
q=input(’Do you desire the muliphase decomposition (Y/N)? ’,’s’); 
if q= езу” | q= шеу” 
пір 
return 
end 


% The Decomposition Routine 


* Determine if we want to ignore the display for each resolution 
* — level's coefficients. 
pltl zinput('Bypass the display of the Coefficients (Y/N)? ’,’s’); 


clc 

phsz2; % This just sets "phs” initially to be NOT equal to 1 ог 0 
phsvct=[]; 

zro- []; 


* The following code sets the number of variables necessary 
% to record the coefficients. ie. lowest=-2, we have 


% c0,cl,c2,d1,d2 
ЕЕ = 

for lvl=-1:-1l:lowest 

while 1 ==1 


phs =input(’Enter the desired Phase [0/1} for this resolution level: ’); 
if phs==1 |phs==0,break,end 
end 
phsvct =[phs phsvct}; 
eval(('c',num2su(abs(lv1)),' 2 phs',num2str(phs), '(c',... 
num2sut(abs(1v1)- 1), ', h);']) 
eval(('d' , num2su(abs(lvl)), ' 2 phs', num2str(phs), '(c',... 
num2str(abs(Ivl)-1),’,g);"}) 


% This if statement runs the display of the coefficients if the flag 
% pltl is set to NO 

if plt] z z'N' | plt] 2 —'n', dspcoef, end 

phs-2; 

eval(( LLL 2 max(LLL,length(c' ,num2sut(-1vl),') *2^ (-1);']) 
zro — [zro 2" (-1v1)-1]; 
end 


% We now do an energy check of the coefficients 


12 


% We now mun the reconstruction option 


reconld 


% END OF WAVID.M 


* Phase Zero decomposition for the 1-D case 

* 

% By: LT J. E. Legaspi 

% Professor Lam 

% US Naval Postgraduate School, Monterey CA 
% e-mail: legaspi@ece.nps.navy.mil 

% lam@ece.nps.navy.mil 

% Phone: (408) 646-3044/2772 


% phsO.m is a function m-file the does a phase 0 decomposition of 
% of the data array using either the h or g coefficients. The main 
% program that calis this routine is "wavld.m" 

% 

% The input arguments for this routine: 


% 

% data -- input data 

% Һ -- the h or g coefficients 

% E -- the length of the data vector 

% H -- the length of the h or g coefficient vector 
ИИИ ае. а ...--.-.--- 
function x 2 phsO(data,h) 

L=length(data); 

H =length(h); 


% The following 4 lines check if the data array is even, for the 
% phase 0 case, we must pad it with one zero if true 
if rem(L/2,floor(L/2)) = = 
data — [data 0]; 
L=L+1; 
end 
data — [zeros(1,H -1) data zeros(1,H-2)]; 
LzL-2*H-3; 
а=0: 
for ii21:2:L-(H-1) 
а=а+1; 
x(1,8) 2 data(ii:ii - H-1)*h'; 
end 
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% Phase One decomposition for the 1-D case 

% 

% Ву: LT J. E. Legaspi 

% Professor Lam 

* US Naval Postgraduate School, Monterey CA 

* e-mail: legaspi@ece.nps.navy.mil 

% lam ece.nps.navy.mil 

* Phone: (408) 646-3044/2772 

Жа... ҚДС eee TRE Ore ere e ur AR E 

* phsl.m is a function m-file the does a phase 1 decomposition of 
% of the data array using either the h or g coefficients. The main 
% program that calls this routine is “wavld.m”" 

% 

% The input arguments for this routine: 


H -- the length of the h or g coefficient vector 


% 

% data -- input data 

* h - the h or g coefficients 

% L -- the length of the data vector 
% 

% 


---ее-----------етее-еетее-е-енкенетет-ете---еенен--ете-ез-ете-н--ен----е-ен--т-- 


function x 2 phsl(data,h) 


L=length(data); 
H =length(h); 


% The following four lines checks if the data vector is odd, if it 
% tme for the phase 1 case, we must pad it with one zero 
if rem(L/2,floor(L/2)) ~ = 0 
data =[data 0); 
L=L+1; 
end 


data — [zeros(1, H-2) data zeros(1,H-2)]; 
L-L-42*H-4; 


а-0; 

for i= 1:2:L-(H-1) 

а=а+1; 

x(1,a) =data(1,ii:ii+H-1)*h’; 
end 


% dspcoef.m Display the coefficients 

% 

% By: LT J. E. Legaspi 

Professor Lam 

US Naval Postgraduate School, Monterey CA 

e-mail: legaspi@ece.nps.navy.mil 
lam@ece.nps.navy.mil 

Phone: (408) 646-3044/2772 


SE OE EE ESSE SS SSS ESSE SSS SET SS SS TS ST SSS SS SST SE ST SEE в б фе чт те би челе ч чечет 


эЧ R R эч эЧ OR 


% dspcoef.m displays the coefficients for each particular 

% resolution level with the appropriate scaling factor. There is 

% also a special routine embedded in this algorithm for the 

* proper bar display of the HAAR case. The main program is "wavld.m" 
% 

* The variables are defined as follows: 

% 

% 1 -- the number of zeros to be padded btwn coeff's 

% 11 -- variable from "w.m" for resoluuon level 
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% fsamp -- sampling frequency of the data 


% wwavelet -- string for the selected wavelet 
С-ы... 
сір 

hold off 


axis('normal') 


% For the HAAR case do the following: 
if wwavelet(1) 2 2 'H' | wwavelet(1) 2 2 'h' 
һ-27(-1м1); 
eval(('ctemp 22" (1vI/2)*c ' , num2su(-Ivl),";']) 
eval(('dtemp —2^" (1IvI/2) *d' .num2str(-0v1),*;']) 
ymax =max(max(ctemp),max(dtemp)); 
if ymax <0,ymax =0;end 
ymin 7 min(min(ctemp), min(dtemp)); 
if ymin > 0,ymin=0;end 
L=length(ctemp); 
if L== 

axis({0 (n*1.5-1)/fsamp 1.2*ymin 1.2*ymax]); 

subplot(211) 

ріоЧ(0 0 n*1.5-1 n*1.5-1]/fsamp,(O ctemp ctemp 0]) 

utle(('Approximation at Resolution level ', num?2strílvl)]) 

xlabel(’Haar wavelet’) 
axis([O (n* 1.5- 1)/fsamp 1.2*ymin 1.2*ymax]); 
subplot(212) 
plot([0 O0 n*1.5-1 n*1.5-1]/fsamp,[Odtemp dtemp 0)) 
üue([' Detail Phase-’ ,num2str(phs)]}) 
pause 


strplt 


else 
x=(0:n:n*L-1]+.5*n; 
x= x/fsamp; 
ахіз((0 x(length(x)) 1.2*ymin 1.2*ymax]); 
subplot(211), bar(x,ctemp) 
uue({’Approximation at Resolution level ', num2str(lv1)]) 
xlabel(’ Haar Wavelet’) 
axis([0 x(length(x)) 1.2*ymin 1.2*ymax]); 
subplot(212), bar(x,dtemp) 
üue([' Detail Phase-' , num2str(phs)]) 
pause 
strpit 

end 


* For NON-HAAR cases, do the following: 
else 
zz2^(-Iv)-1; 
eval(('ctemp 22" (1Ivl/2)*ptzero(c' ,num2su(-10v1),' .z);']) 
eval(('dtemp 22^(1vI/2)* ptzero(d' num2sti(-1vl),' ,z);']) 
ymax — max(max(ctemp),max(dtemp)); 
if ymax < 0,ymax =0; end 
ymin = min(min(ctemp),min(dtemp)); 
if ymin > 0,ymin=0;end 


% We can display the data with a window size of the original 
% input, or show the spillover values 
q=[{ Enter the type of display:’ 


' 1 - Windowed Display ' 
' 2- Non-Windowed ']; 
сіс 
disp(q) 
disp(  ']) 
q=input(’ ’); 
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if q= =2 

x=[0:length(ctemp)-1]+.5; 

x= x/faamp; 

axis([O x(length(x)) ymin ymax]); 

subplot(211), bar(x,ctemp) 

uile([’ Approximation at Resolution level ',num2astr(lvl)]) 

xlabel([wwavelet,’ Wavelet’]) 

ах18([0 x(length(x)) ymin ymax]); 

subplot(212),bar(x,dtemp) 

uue([' Detail Phase- ', num2su(phs)]) 

xlabel(['Non-Windowed display -- max windowed value: ’,... 
num2str{length(x)- 1)]) 

pause 

strplt 

else 

x=[0:LL-1] +.5; 

x=x/fsamp; 

axis({0 x(length(x)) ymin ymax]); 

subplot(211),bar(x,ctemp(1,1:LL)) 

tide([' Approximation at Resolution level ',num2su*(lvl)]) 

xlabel([wwavelet," Wavelet']) 

axis({0 x(length(x)) ymin ymax]); 

subplot(212),bar(x,dtemp(1,1:LL)) 


ше(Г Пеш Phase- ',num2str(phs)]) 
xlabel(['Windowed display -- max windowed value: ',num2str(LL-1)]) 
pause 
strplt 
end 


end 
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% enrgid.m Determining the Energy in cach resolution level 


% One Dimensional Case 

% 

% By: LT J. E. Legaspi 

% Professor Lam 

% US Naval Postgraduate School, Monterey CA 
% e-mail: legaspi@ece.nps.navy.mil 

% lam@ece.nps.navy.mil 

% Phone: (408) 646-3044/2772 


ne акта ата т а та Иа ааа 


% For the single phase case, the energy is normalized to resolution 
% level 0. All the energies of the "c'a" and the "d's" of a 

% particular level should add up to the energy of the "c's" of the 
% next higher level. 


enrgc -- energy array of the "c's" 

enrgd -- energy array of the "d'a" 

norm -- energy of the data "cO" 

esum -- energy sum array (enrgc-t enrgd) 


a& a& aS 33 эЧ эЧ 


* The calling routine i$ "wavld.m" . 
7; atrplt.m is a necessary m-file 


enrgc —zeros(1,-lowest); 
enrgd —zeros(1 ,-lowest); 
а-0; 
norm —sum(c0 .^2); 
for ii=lowest:1:-1 
а=а+1; 
eval(['enrgc(a) 2sum(c' ,num2str(-ii), ^2); ']) 
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eval(['enrgd(a) 2 sum(d',num?2su(-i),'.^2);']) 
end 
enrgc — enrgc/norm; 
епгрс = [enrgc 1]; 
enrgd 2 enrgd/norm; 
clg 
axis([lowest-2 1 0 1.2)); 
subplot(211),bar([lowest:0), [enrgc]) 
title({’ Normalized energy of the "c" coefficients Phase: ’,... 
setstr(flipIr(phsvct + 48))]) 
xlabel(’ Resolution level’) 
axis([lowest-2 1 0 1.2]); 
subplot(212),bar([lowest:-1],[enrgd]) 
utle(’ Normalized energy of the "d" coefficients) 
xlabel(' Resolution level") 
pause 
strpit 


axis 
% Energy check of the previous plot 


esum = епгрс + [enrgd 0]; 


clc 

dd 2[' Venfication of the energy in each resolution level ' 
"Level Energy in Energy in C Difference’ 
: level m-1 level m , 
i c+d е 
К Тр 

disp(dd) 

j=-lowest; 


for i= 1:-lowest 
а = езшт(1);Ь = епгес(1,1+1);с = 1-j;d =abs(a-b); 
fprintí(' 2.07 %17.15/ ЖЫНДАР "сар 
fprintf(’ %17.15f\n’,d) 
]=]-1; 

спа 

раиве 

сіс 

axis; 


% reconid.m One Dimensional Single Phase Reconsruction 

% 

* By: LT J E. Legaspi 

Professor Lam 

US Naval Postgraduate School, Monterey CA 

e-mail: legaspi@ece.nps.navy.mil 
lam@ece.nps.navy.mil 

Phone: (408) 646-3044/2772 


3% 38 38 3% 28 3% OR 


% This routine conducts a reconstruction of the the data with the 
% proper phase selected between resolution levels. This is a sub- 
% routine for the main wavelet program "wavld.m". 

% 

% The following variables are further defined: 


% 

%  hh& gg -- reconstruction coefficients 

% lowest -- defined in "wavld.m" 

Ф смотк -- reconsuucted c'values at a particular lvl 
% ^ phsvct -- phase vector defined in "wavld.m". 


TI 


The following m-files are necessary for proper operation: 


* 

% 

% 

% wavld.m -- main program 

% rphs]l.m -- reconstructon phase- 1 function m-file 
% rphsÜ0.m  -- reconstruction phase-0 function m-file 
*  strplt.m -- storage of the displays 

% 
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cle; 

disp([’ Т? 

q = input(’Do you desire to do the Reconstrucuon Comparison (Y/N)?’,’s’); 
cle 


fq ИУ. 

% generation of the reconstruction coefficients 

hh = fliplr(b);gg —fliplr(g); 

i=]; 

eval([('cwork —c',num2str(-lowest),';']) 

for lv1=lowest:1:-1 
eval(['dwork 2d',num2st(-1v)),';'] 
eval(['cwork rphs' ,num2su(phsvct(i)), '(cwork,dwork,hh,gg);']) 
eval(['a -length(c' ,num2astr(-1vl- 1), 5;']) 
cwork — cwork(1:8); 
u=utl; 


phswrk 7 phsvct(1:length(phsvct)-1); 

* Comparison of the Original and the Reconstructed data 

clg 

subplot(211) 
eval(['bar(c',num2str(-1vl-1),')' D 
utle(['Original Data Phase: ',sestr((phswrk - 48))p) 
xlabel([' Wavelet: ',wwavelet]) 

subplot(212) 

bar(cwork); 

utle(' Reconstructed Data") 
xlabel(('Resoluuon level ',num2str(lvl-- 1)]) 

pause 

strplt 

clg 

a —- eval(['c' ,num2str(-1vl- 1), ' -cwork' p; 

bar(abs(a)); 

utle(' Absolute Error in the Reconstruction") 

xlabel([' Resolution Level ',num2str(lvl *- 1)]) 

pause 

suplt 

clg 

phswrk z phswrk(1:length(phswrk)-1); 

end 
end 


function out =rphs0(cwork,dwork,bh, gg) 


Ез КЫ ӨЕ КЫЛЫК ЕКИ КУКУ» ККК КУК 


?; rphs0.m Reconstruction Phase Zero routine 


% 
% By: LT J. E. Legaspi 
% Professor Lam 


% US Naval Postgraduate School, Monterey CA 
% e-mail: legaspi@ece.nps.navy.mil 

% lam@ece.nps.navy.mil 

% Phone: (408) 646-3044/2772 

% 


*&$*5$9$$9***$$$*$***x***x*****$$$$$$$$*$$$*******$$$*$t*$t*tt$*$t*t********$t** 
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% 

% This routine does a reconstruction of a particular level with a 
% phase 0 selection. The calling program is "*wavld.m". Input 
% arguments are: 


ж 


% cwork -- *c* coefficients of a particular level 
% dwork -- 'd* coefficients of a particular level 
% hh -- reconstruction coeff’s for the °c’s” 

% gg  -- reconstruciton coeff's for the "d's* 
% 
% 
% 


"out" is the output argumemt that equals the c coeff's of the 
next higher resolution level. 


L=length(cwork); 
хі-П; 
х2 =[; 
for j=1:L 
x1 =[x1 cwork(1,}) 0]; 
x2 z[x2 dwork(1,j) 0]; 
end 
x1 2 [x1 zeros(1,length(hh)-2)]; 
x2 2 [x2 zeros(1, length(gg)-2)]; 


а=0; 
for j 2 1:length(x1)-(length(hh)-1) 
а=а+1; 


out(] ,a) 2 x1(1,j:j * length(hh)- 1) *hh' *- x2(1,j:j t length(gg)-1)*gg '; 
end 


function out 2 rphsl(cwork dwork hh, gg) 


К © © ооо ооо ооо оо ооо оо ооо ооо зо 


% rphsl.m Reconstruction Phase One routine 


% 

% By: LT J. E. Legaspi 

% Professor Lam 

% US Naval Postgraduate School, Monterey CA 
% e-mail: legaspi@ece.nps.navy.mil 

% lam@ece.nps.navy.mil 

% Phone: (408) 646-3044/2772 


ооо ооо ооо сос I HOLS ооо ооо соо 


% 

% This routine does a reconstruction of a particular level with a 
% phase 1 selection. The calling program is “wavld.m". Input 
% arguments are: 


cwork -- “°c” coefficients of a particular level 
dwork -- 'd' coefficients of a particular level 
hh  -- reconstruction cocff's for the "c's" 
£g  -- reconstruciton coeff's for the "d's" 


“out” is the output argumemt that equals the c coeff's of the 
next higher resolution level. 


азая Яя GR OS OR 


L=length(cwork); 
х1=[]; 
x2 z []; 
for ;=1:L 
x1 =[x1 cwork(1,}) 0); 
x2 =[x2 dwork(1,j) 0); 
end 
xl =[0 x1 zeros(1,length(hh)-2)]; 
x2 — [0 x2 zeros(1,length(gg)-2)]; 
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а=0; 
for j 2 1:length(x1)-(length(hh)- 1) 

а=а+]; 

out(1,2) 2 x1(1,j:j * length(hh)- 1)*hh' * x2(1,j:j t length(gg)-1)*gg'; 
end 


q*ss45552295999992999€89989€99999999999»9 99929 92999999 99999 2» 22 99.9 


% plane.m develop the phase plane diagram 


% 

% By: LT J. E. Legaspi 

% Professor Lam 

% US Naval Postgraduate School, Monterey CA 
* e-mail: legaspi@ece.nps.navy.mil 

% lam@ece.nps.navy.mil 

% Phone: (408) 646-3044/2772 


pA S LL LS e LaL LO Sa i d d od ad ao I Lo s ee eee аа 


% Phase Plane Diagram for the One Dimensional Single Phase case 
% This routine is called by the parent routine "wavld.m^. All the 
% variables defined in "wavld.m” are used here. New variables: 

% 

% c --(l-lowest) X LL "c" coefficient matrix 

% d --(lowes) X LL "d" coefficient matrix 

%  N,NI -- variables for the number of contour levels 

% 

% The main subroutine used is ptzero.m funtion m-file 


cle 
clg 


% Setup the c and d matrices from the "c*" and "d*" variables 


c -zeros(-lowest 4 1,11); 

d =zeros(-lowest, LL); 

c(1,1:length(c0)) 2 c0; 

clear cwork dwork 

for lvl=1:-lowest 
cwork =eval((’ptzero(c’ ,num2str(Iv]),’ ,zro(Ivl))’}); 
dwork —eval(('ptzero(d'  num2str(lvl), ', zro(1vD))']); 
c(Iv1 1,1:length(cwork)) 2 cwork; 
d(Ivl, 1: length(dwork)) =dwork; 


end 

c=flipud(c.*2); 

d 7 flipud(d. ^2); 

[i j] 7 size(c); 

х=0:)-1;ус=0:1-1;уа=1:1-1; 

mesh(c); 

title([’Energy Display of the Approximation Phase: ',... 
setsu(fliplr(phsvct 4- 48))]); 

pause 

N=10; 


contour(c, 10,x,yc); 
utle(['Contour Display of the Approximation Phase: ’,.. 
setsu(fliplr(phsvct 4- 48))]); 
ylabel(’-n resolution level’); 
xlabel(’10 contour levels’) 
pause 
clg 
ql = input(’ Do you desire to zoom in on the display (Y/N)?’,’s’); 
while ql=='Y’ | ql == 'y', 
q2 7 input(' Do you want to see thc plots again (Y/N)?’,’s’); 
if 42- =’y’ | q2 == ty”. 
mesh(c);title(' Previous Energy Display of the Approximation’); 
pause 
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contour(c, N , x, yc);title(' PreviousContour Display of the Approximation’); 
ylabel(’-n resolution level’); 
xlabel(['Number of countours: ',num2str(N)]) 
pause 
end 
clc 
вр Т) 
disp(['x range: 0 to ',num?2str(j-1)]) 
disp(['y range: 0 to -',num2str(i-1)];disp([. — ]) 
x] =input(’Enter the minimum x-value:’)+ 1; 
x2 =input(’Enter the maxumum x-value:’) +1; 
y]=input(’Enter the minimum “magnitude” y-value (least negative):’)+1, 
y2 = input( Enter the maxumum "magnitude" y-value (most negative):') * 1; 
yl -abs(yl);y2 —-abs(y2); 
у= 1:1:;у = ріку); 
у11 =у(у2); 
у22 =у(у1); 
сір 
mesh(c(y11:y22,x1:x2)) 
ütle(['Zoomed Energy Display of the Approximation Phase: ’,... 
setsu(fliplr(phasvct 4- 48))]); 
pause 
strplt 
clg 


N1=10; 
while N1 < 999 
М-М1; 
contour(c(y11:y22,x1:x2), N1,x1:x2,y1-1:y2-1) 
title(['Zoomed Contour Display of the Approximation Phase: ’,... 
setsu(fliplr(phsvct4-48))]); 
ylabel(’-n resolution level’) 
xlabel(['Number of contours: ',num2su(N1)]) 
pause 
вирі 
N1 z input(' Enter the number of contour levels (999 to continue) °); 
end 
ql =input(’Zoom in or out Further (Y/N)? ’,’s’); 
end 
cle 


% Now the Detail Signal 


mesh(d); 

title([’Energy Display of the Detail Phase: ',setsu(fliplr(phsvct - 48))]); 
pause 

clg 


contour(d,10,x,yd); 

title(['Contour Display of the Detail Phase: ',setstr(fliplr(phsvct4-48))]); 
ylabel(’-n resolution level’); 

xlabel(' 10 contour levels’) 


pause 
clg 

clc 

91 z input('Do you desire to zoom in on the display (Y/N)?’,’s’); 
while q12 Z'Y'| ql == у’, 


q2 — input(' Do you want to see the plots again (Y/N)?','8*); 
if q222'Y'|q2 2-'y', 
mesh(d);ttle('Previous Energy Display of the Detail’); 
pause 
contour(d, N, x, yd);title(' PreviousContour Display of the Detail’); 
ylabel(’-n resolution level’); 
xlabel(['Number of contours: ',num2su(N)]) 
pause 
end 


8l 


disp(['x range: 0 to ',num2su(j-1)]) 
disp(['y range: -1 to -',num2su(i-1)]);disp([. ]) 
xl 2 input('Enter the minimum x-value:") * 1; 
x2 =input(’Enter the maxumum x-value:’)+1; 
y1 =input(’ Enter the minimum “magnitude” y-value (least negauve):"); 
if yl == 0, yl=1; end; 
y2 =input(’Enter the maxumum “magnitude” y-value (most negative): "); 
y 1 =abs(y1);y2 =abs(y2); 
у=1:1:1-1;у=ЯрЩуУ); 
у11-у(у2); 
у22=у(у1); 
сір 
mesh(d(y11:y22,x1:x2)) 
utle(['Zoomed Energy Display of the Detail Phase: ’,... 
setstr(flip]r(phsavct+48))]); 
pause 
strplt 
сір 
while N « 999 
contour(d(y11:y22,x1:x2), N,x1:x2,y1:y2) 
title([’ Zoomed Contour Display of the Detail Phase: ’,... 
setstr(fliplr(phsvct *- 48))]); 
ylabel('-n resolution level") 
xlabel(['Number of contours: ',num2su(N)]) 
pause 
strplt 
= input’ Enter the number of contour levels (999 to end): "); 
end 
q1 2 input(' Zoom in or out further (Y/N)? ’,’s’); 
end 


% Multiphase routine for The Discrete Wavelet Decomposition Routine 
% 

% By: LT J. E. Legaspi 

% Professor Lam 

% US Naval Postgraduate School, Monterey CA 

% e-mail: legaspi@ece.nps.navy.mil 

% lam@ece.nps.navy.mil 

% Phone: (408) 646-3044/2772 

% 


This routine does a multiphase discrete wavelet decomposition of 
a one-dimensional data vector. This routine requires the following 
M-files: 
wavld.m 
strplt. m 
zÜpad.m 
maxSamp —(Nh-1)*(2* (-lowest)-1) - LL; 
coeffc 2 zeros(-lowest -- 1], mmaxSamp); 
coeffc(-lowest - 1,1: LL) 2 c0; 
coeffd —zeros(-lowest,mmaxSamp); 
numcoef —zeros(-lowest * 1,1); 


a& a& aq a a8 33 aX aq 


Ф *** The Main Multi-Phase Decomposition Algorithm *** 
numcoef(-loweat * 1) Z LL; 
11; 
for row =-loweat:-1:1 
numcoef{row) =(Nh-1)*(27i-1)+ LL; 
pl =zOpad(coeffc(row + 1,1: numcoef(row + 1)),h,i); 
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coeffc(row,1:length(p1)) 2pl; 
p2 7zOpad(coeffc(row + 1,1:numcoef(row * 1)),g.1); 
coeffd(row,1:length(p2)) 2 p2; 
clear pl p2 
1=1+1; 
сікі 


% Display of Coefficients at each resolution level 
clc 
q=input(’Do you desire to see values at diff resolution levels (Y/N)? ’,’s’); 
if q — — dd | q == у’ 
plot([0:LL-1],c0,' **) 
xlabel(['n sampling points n*' num2str(1/fsamp),' seconds']); 
ylabel(’ Amplitude’) 
uUe(’Original Data (resolution level 0) Multiphase case’); 
pause 
strplt 
clg 
index 7 [-lowest:- 1:1]; 
x1=0;x2=-lowest+1,; 
cle 
disp([’ Levels to pick are from -1 to ',num2su(lowest)]) 
while x1« 1 | x2» -lowest 
x1 =input(’Enter first resolution level of the desired range: ”); 
x2 =input(’Enter second resolution level of the desired range: ’); 
x1 =abs(x1);x2 =abs(x2); 
if x1 > x2, temp=x2; x2=x1; xl =temp,; end 
end 
while x1 < =x2 
pts =numcoef(index(x1)); 
Tc- coeffc(index(x1), 1:pts); 
Td =coeffd(index(x1),1:pts); 
minpts = min(min(Tc),min(Td));ifminpts > 0,minpts = 0; end 
maxpts = max(max(Tc),max(Td));ifmaxpts < 0,maxpts =0;end 
if minpts = =0 & maxpts= =0, maxpts =.5;end 
u=(0 pts 1.2*minpts 1.2*maxpts}; 
axis(u); 
subplot(211).plot([0:pts-1], coeffc(index(x1), 1:pts),'*") 
title([’ Multiphase Approximation Coeff at Resolution Level -' , num2su(x1)]) 
axis(u); 
subplot(212), plot([0:pts-1], coeffd(index(x1), 1:pts),'*") 
title(’ Multiphase Detail Coefficients’) 


xlabel([’n sampling points n*’ num 2str(1/fsamp),’ seconds’); 
pause 
вирі 
clg;clc 
axis; 
х1=х1+1; 
епа 
епа 
лс 
% Multi-Phase Energy Check 
enrgldmp 
ЕЕ аеоси еее онооно 
% Phase Plane Determination 


c=coeffc.*2; 
d=coeffd.*2; 

[1 j] 2 size(coeffc); 
x=0:LL-1; 

ус=0:1-1; 

yd=1:1-1; 

cc 2c(1:-lowest - 1,1: LL); 
clc 

mesh(cc); 

title(’ Multiphase Energy Display of the Approximation’) 
pause 


suplt 
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clg 
N=10; 
contour(cc, 10,x, yc); 
uile(’ Multiphase Contour Display of the Approximation’); 
ylabel(’-n resolution level’); 
xlabel(’10 contour levels’) 
pause 
strplt 
сір 
41 7 input('Do you desire to zoom in/out on the display (Y/N)?’.’s’), 
while ql==’Y’ | ql == 'y', 
q2=input(’Do you want to see the plots again (Y/N)?’,’s’); 
if g2==’Y’ | q2 =='y’, 
mesh(cc) 
uue(’Previous Multiphase Energy Display of the Approximation’), 
pause 
contour(cc,N,x,yc) 
title(" Previous Multiphase Contour Display of the Approximation’); 
ylabel(’-n resolution level’), 
xlabel({’ Number of contours: ’,num2str(N)}) 
pause 
end 
cle 
disp ']’) 
disp(['x range: 0 to ',num2str(j-1)]) 
disp(['y range: O to -',num2sti(i-1));disp( СҮ) 
x] =input(’ Enter the minimum x-value:’) +1, 
x2 = input(’ Enter the maxumum x-value:") + 1; 
yl =input(’Enterthe minimum "magnitude" y-value (least negative):') +1; 
y2=input(’Enter the maxumum “magnitude” y-value (most negative):’) +1; 
y 1 =abs(y1);y2 =abs(y2); 
yzl:l:uy-flplr(y); 
у11-у(у2); 
у22-у(у1); 
clg 
сс-с(у11:у22,х1:х2); 
mesh(cc) 
utle(’Zoomed Multiphase Energy Display of the Approximation’); 
pause 
strplt 
clg 
N=10: 
М1 =10; 
while N1« 999 
NzNI; 
contour(cc, N1,x1:x2,y1-1:y2-1) 
utle(’Zoomed Multiphase Contour Display of the Approximation’), 
ylabel(’-n resolution level’) 
xlabel({’ Number of contours: ',num2str(N1)]) 
pause 
strplt 
N1=input(’ Enter the number of contour levels ( & RET? to continue) ’); 
end 
ql =input(’Zoom in or out Further (Y/N)? ’,'s’); 
end 
cle 


% Now the Detail Signal 

dd =d(1:-lowest,1:LL); 

mesh(dd) 

utle(’ Multiphase Energy Display of the Detail’) 
pause 

strplt 

clg 


contour(dd, 10,x,yd) 
utle(' Mulüphase Contour Display of the Detail); 
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ylabel(’-n resolution level’); 
xlabel(’10 contour levels’) 
pause 

strplt 

clg 


clc 
ql =input(’ Do you desire to zoom in/out on the display (Y/N)?’,’s’); 
while ql=='"Y" | ql == у’, 
q2 2 input('Do you want to see the plots again (Y/N)?’,’s’); 
if q= =Y’ | 42 =='у', 
mesh(dd);utle(' Previous Multiphase Energy Display of the Detail); 
pause 
contour(dd, N,x,yd);title(’ PreviousMultiphase Contour Display of the Detail’), 
ylabel("-n resolution level"); 
xlabel([' Number of contours: ',num2str(N)]) 
pause 
end 
disp(['x range: 0 to ',num2str(j-1)]) 
disp(['y range: -1 to -",num2str(i-1)));disp([ — ]^) 
x1=input(’ Enter the minimum x-value:') * 1; 
x2 zinput('Enter the maxumum x-value: ') * 1; 
ylzinput('Enter the minimum "magnitude" y-value (least negative): '); 
if yl == 0, yl=1; end; 
y2=input(’Enter the maxumum “magnitude” y-value (most negative): "); 
yl-abs(yl);y2-abs(y2); 
у=1:1:1-1;у =ЯфрЩу); 
У11-у(у2); 
у22 =у(у1); 
сір 
dd =d(y11:y22,x1:x2); 
mesh(dd) 
utle(’Zoomed Multiphase Energy Display of the Detail’); 
pause 
surplt 
clg 
while N « 999 
contour(dd,N,x1:x2,y1:y2) 
title('Zoomed Mulüphase Contour Display of the Detail’); 
ylabel('-n resolution level’) 
xlabel([' Number of contours: ',num2str(N)]) 
pause 
surplt 
=input(’Enter the number of contour levels (< RET > to end): '); 
end 
q1 zinput('Zoomin or out further (Y/N)? ’,’s’); 
end 


mprecon 
Е 12га 
% The End 
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% ежежететевееееееееетеететезтеетекететететкеттен етене тет al 


% enrgldmp.m Determining the Energy in each resolution level 


% Multi-Phase One Dimensional Case 
% 

* By: LT J. E. Legaspi 

% Professor Lam 

% US Naval Postgraduate School, Monterey CA 

% e-mail: legaspi@ece.nps.navy.mil 

% lam@ece.nps.navy.mil 

% Phone: (408) 646-3044/2772 


Жыннан 


% For the Multi-phase case, the energy is normalized to resolution 
% level 0. All the energies of the "c's" and the 'd's" ofa 

% particular level should add up to the energy of the "c's" of the 

% next higher level. The factor of two is accounted for in the Multi- 
% Phase case 

enrgc -- energy array of the “с?” 

enrgd -- energy array of the 74787 

norm -- energy of the data "cO" 

esum -- energy sum array (enrgc-t enrgd) 


R 8 a3 a& aq av 


% The calling routine is "mp.m" . 
% wavid.m and strplt.m are also necessary m-files 


епгес-?егов(1,-Іоу/ез1); 
enrgd =zeroa(1 ,-lowest); 
a=0; 
norm=sum(c0 .*2); 
for row =-lowest:-1:1 
а=а-+ 1; 
enrgc(row) 2 sum(coeffc(row,:).^2)*2^(-8); 
enrgd(row) 2 sum(coeffd(row,:).^2)*2^ (-8); 
end 
enrgc - enrgc/norm; 
enrgc 7 [enrgc 1]; 
enrgd 7 enrgd/norm; 
clg 
axis([lowest-2 ] 0 1.2)); 
subplot(21 1),bar([lowest:0]. [enrgc]) 
utle([’ Normalized energy of the "c" coefficients Multi-Phase Case’]) 
xlabel(’ Resolution level’) 
axis([lowest-2 1 0 1.2]); 
subplot(212), bar([lowest:-1), [enrgd]) 
title(’Normalized energy of the "d" coefficients") 
xlabel(’ Resolution level’) 
pause 
strplt 


axis 


% Energy check of the previous plot 


esum =enrgc + [enrgd 0); 
cle 
dd=[’ Verification of the energy in each resolution level — ' 
"Level Energy in Energy in C Difference’ 
level m-1 level m ” 
c+d f 
p ull le e LS " 
disp(dd) 
j=-lowest; 


for i=]:-lowest 
a =csum(i);b=enrgc(1,i+ 1);c=1-j;d =abs(a-b); 
fprintí(' 992.0f — 9617.15f 9?017.15f ’,c,a,b) 
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fprintf(’%17.15f\n’,d) 
ізі; 

епа 

pause 

clc 

clg 

axis; 


кт ккккж т же жетке итте жекке ееткке ә өөежжеееетже 


% mpplane.m develop the phase plane diagram 


% 

% Ву: LT J. E. Legaspi 

% Professor Lam 

% US Naval Postgraduate School, Monterey CA 
% e-mail: legaspi@ece.nps.navy.mil 

% lam@ece.nps.navy.mil 

% Phone: (408) 646-3044/2772 


Ж 7777758 %9996%%566606Фее6600000666600000000006060ө000060ө6000406 


% Phase Plane Diagram for the One Dimensional Single Phase case 
% This routine is called by the parent routine “wav}d.m*. All the 
% variables defined in “wavld.m* are used here. New variables: 

% 

*$ c --(l-lowest) X LL "c° coefficient matrix 

*$ d --(-lowest) X LL "d" coefficient matrix 

%  N,NI - variables for the number of contour levels 

% 

% Тһе main subroutine used is ptzero.m funtion m-file 


cle 
clg 


% Setup the c and d matrices from the "c** and "d*" variables 


c zzeros(-lowest- 1, LL); 

d zzeros(-lowest, LI); 

c(1,1:length(c0)) 2 c0; 

clear cwork dwork 

for lvl=1:-lowest 
cwork — eval(['ptzero(rmpc' ,num2str(lv1),', 2^ (1v])-1)]); 
dwork - eval(('ptzero(rmpd' ,num2st(1v1),',2^ (1vD-1)']); 
c(Iv14 1,1:length(cwork)) 2 cwork; 
d(lvl,1:length(dwork)) 2 dwork; 

end 

c z flipud(c.^2); 

d — flipud(d. ^2); 


[i 3] =size{c); 
х=0:}-1;ус=0:1-1;у9=1:1-1; 
mesh(c); 
title(('* MP* Energy Display of the Approximation Phase: ’,... 
setsu(fliplr(phsvct - 48))]); 
pause 
N z10; 
contour(c,10,x,yc); 
tile(( * MP* Contour Display of the Approximation Phase: ’,... 
setstr( fliplr(phsvct + 48))]); 
ylabel(’-n resolution level’); 
xlabel(' 10 contour levels’) 
pause 
clg 
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q1 zinput(' Do you deaire to zoom in on the display (Y/N)?’,’s’); 
while q1 2 2'Y' | 91 == 'у’, 
q2 7 input(' Do you want to see the plots again (Y/N)?’,’s’); 
if q22 2'Y' | q2 == Vin 
mesh(c); 
title(’*MP* Previous Energy Display of the Approximation’); 
pause 
contour(c,N,x,yc); 
uue(’*MP* Previous Contour Display of the Approximation’); 
ylabel(’-n resolution level’); 
xlabel({’ Number of countours: ',num2str(N)]) 
pause 
end 
cle 
disp( УТ) 
disp(['x range: 0 to ’ ,num2str{j-1)]) 
disp(['y range: O to -",num2str(i-1));disp([ — ]?) 
x1 =input(’ Enter the minimum x-value:') + 1; 
x2 =input(’Enter the maxumum x-value:’) +1; 
y 1 =input(’ Enter the minimum “magnitude” y-value (least negative):’) +1, 
y2=input(’Enter the maxumum "magnitude" y-value (most negative):) * 1; 
y 1 =abs(y1);y2 =abs(y2); 
у=1:1:;у = ріку); 
у11-у(у2); 
у22 =у(у1); 
сір 
mesh(c(y11:y22,x1:x2)) 
ütle([' *MP* Zoomed Energy Display of the Approximation Phase: ’,... 
setstr(fliplr(phsvct 4-48))]); 
pause 
suplt 
clg 


N1=10; 
while N1« 999 
NzNl; 
contour(c(y1ll:y22,x1:x2), N1,x1:x2,y1-1:y2-1) 
ütle([' * MP* Zoomed Contour Display of the Approximation Phase: ’,... 
setstr(fliplr(phsvct+48))]); 
ylabel(’-n resolution level’) 
xlabel({’ Number of contours: ’,num2str(N 1)]) 
pause 
strplt 
N1 =input(’ Enter the number of contour levels (999 to continue) ”); 
end 
ql =imput(’Zoom in or out Further (Y/N)? ’,’s’); 
end 
cle 


% Now the Detail Signal 


mesh(d); 

utle([ * MP* Energy Display of the Detail Phase: ’,... 
setsu(fliplr(phsvct4- 48))]); 

pause 

clg 


contour(d, 10, x, yd); 

utle([ *MP* Contour Display of the Detail Phase:',... 
setsu(fliplr(phsvct 4- 48))]); 

ylabel(’-n resolution level’); 

xlabel(’10 contour levels’) 

pause 

clg 


cle 
q1- input('Do you desire to zoom in on the display (Y/N)?’,’s’); 
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while ql==’Y’ | ql == уу”, 
q2 =input(’Do you want to see the plots again (Y/N)?','s'); 
if q2==’Y’ | q2 = ='у’, 
mesh(d);title('Previous Energy Display of the Detail’); 
pause 
contour(d, N,x,yd);title(' PreviousContour Display of the Detail’); 
ylabel(’-n resolution level’); 
xlabel([' Number of contours: ',num2str(N)]) 
pause 
end 
disp(['x range: 0 to ',num2str(j-1)]) 
disp(['y range: -1 to -',num2sut(i-1)]);disp([ ]’) 
x1 zZinput('Enter the minimum x-value:') * 1; 
x2 =input(’Enter the maxumum x-value:’) +1; 
y 1 =input(’Enter the minimum "magnitude" y-value (least negative):'); 
if yl == 0, yl=1; end; 
y2=input(’Enter the maxumum “magnitude” y-value (most negative):’); 
y1=abs(y1);y2 =abs(y2); 
y =1:1:4-1;,y =flipir(y); 
yll=y(y2); 
у22 =у(у1); 
сір 
mesh(d(y11:y22,x1:x2)) 
ütle([ * MP* Zoomed Energy Display of the Detail Phase: ',... 
setstr(fliplr(phsvct 4- 48))]); 
pause 
вирі! 
сір 
while N « 999 
contour(d(y11:y22,x1:x2), N,x1:x2,y1:y2) 
title(( * MP* Zoomed Contour Display of the Detail Phase: ’,... 
setsu(fliplr(phsvct 4- 48))]); 
ylabel(’-n resolution level’) 
xlabel([(' Number of contours: ',num2su(N)]) 
pause 
strplt 
N =input(’Enter the number of contour levels (999 to end): ’); 
end 
ql =input(’Zoom in or out further (Y/N)? ’,’s’); 
end 


function x =zOpad(Data,h, reslvl) 

ee 2 2 SPSS OES AS SEES LSE SATOH SEA ESE SSSR SSSA SE SS SESSESS SOLES SESS 
% zOpad.m — *** For the Multi-phase decomposition routine *** 

% One Dimensional Case 

% 

% By: LT J. E. Legaspi 

% Professor Lam 

% US Naval Postgraduate School, Monterey CA 

% e-mail: legaspi@ece.nps.navy.mil 

% lam@ece.nps.navy.mil 

% Phone: (408) 646-3044/2772 

% 
% 


Фескеооеоокееоееоеоеевеооевеооеоооеееоееееосеиежкекооеоеоеоеееееоееееегеге 


% This routine makes a row vector of the lower coefficients of 
% a particular resolution level. Use for only the Multi-phase 
% decomposition One-Dimensional case !! 


Data is the data vector 


% 
% 
% h is either the h coeff or the g coeff vector 
% realv] is the resolution level 
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m — reslvl; 

L=length(h); 

ізі 

x=0; 

for i=1:L 
x=x+h(j)*{zeros(1,(i-1)*2*(m-1))Data zeros(1,(j-1)*2*(m-1))); 
j=j-1; 

end 


function binary =num2bin(number, length1) 


ME Ек. лал ОА Оаа ла НЫ 
% NUM2BIN.M Number to Binary conversion 
% 

% By: LT J. E. Legaspi 


% Professor Lam 

% US Naval Postgraduate School, Monterey CA 
% e-mail: legaspi@ece.nps.navy.mil 

* lam@ece.nps.navy.mil 

% Phone: (408) 646-3044/2772 

% 
% 


% This routine conducts a conversion of a base 10 integer to a base 2 
% data array i.e. [10101 1} 

% 

% The following variables are further defined: 


% 
% binary” -- binary output data array 
% length] -- input desired length of the display 
% number -- base 10 integer 
% 
% The following m-files are necessary for proper operation: 
% 
% wavid.m -- main program (great grandparent routine) 
% mpm -- grandparent routine 
* mprecon.m -- multiphase reconstuction calling routine 
Eee em wee EE E EE EE E E EE EE E E RE E E E E E ERR EE EUER Em EE m 
binary —[]; 
while number > .5 

number =number/2; 


if number-fix(number) = =. 
number =number-.5; 
binary -[1 binary]; 


else 
binary =(0 binary}; 
end 
end 
if length(binary) < length] 
binary = [zeros(1,abs(lengthl-length(binary)))binary]; 
end 
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function num =bin2num(bin) 


ENERO S soos beeen oraaa 
% BIN2NUM.M Binary to number conversion 
% 

% Ву: LT J. E. Legaspi 

% Professor Lam 


% US Naval Postgraduate School, Monterey CA 
* e-mail: legaspi@ece.nps.navy.mil 

% lamGece.nps.navy.mil 

% Phone: (408) 646-3044/2772 

% 


% This routine conducts a conversion of a base 2 data array 
% i.e. [1010 1 1) to a base 10 integer 

% 

* The foliowing variables are further defined: 

* 

Ф Ып -- binary output data array 

%  num_ -- base 10 integer 

% 


% The foliowing m-files are necessary for proper operation: 


% wavid.m_ -- main program (great grandparent routine) 
% mpm -- grandparent routine 
%  mprecon.m-- muluphase reconstuction calling routine 


пит = 0; 
8 — length(bin); 
power-a; 
for ii—1:a 
num —num 4 bin(ii) *2^ (power-ii); 
end 


9] 


B. TWO-DIMENSIONAL ROUTINES 


ООРУУ ФУУУ ИЕЫ ЕКЕЖ ЕКЕУ ҮР УКУ КУ КУЕ o "ev T C eT e 
% WAV2D.M 15 SEP 92 

Ф Wavelet decomposition algorithm for the two-dimensional case 

% for either the single selectable (of 4 possible) or the 

* multiple phase case. 


% 

* By: LT J. E. Legaspi 

% Professor Lam 

% US Naval Postgraduate School, Monterey CA 
% e-mail: legaspi@ece.nps.navy.mil 

% lam(ece. nps.navy.mil 

* Phone: (408) 646-3044/2772 


ОЕК ӨК ӨР ИЕЫ ЕЕ УУРУУ РҮ РИ УУУУ УКИ КУУРУ КУЕ ККУУ 


* This routine does a two-dimensional dicrete wavelet decomposition 
* of a two-dimensional data array. The following routines are also 
% necessary for the proper operation of the algorithm: 


chkputer.m -- checks for the minimum hardware 
mp2d.m -- 2-D muluple phase decomposition 
pzro2d.m -- 2-D row & col zero padding 

suplL.m  -- metafiles of selected plots 

hcoefin — -- routine for selected h coefficients 

phas ?? .m -- routines for phase 00,01,10,11 decomp 
enrg2d.m  -- energy check for the 2-D case 
recon2d.m -- 2-D reconstruction 


———  À өз в «в «в өз «е өз «сез өз өз өз «в «г <> өз <> өз «> «> ез ез өз «> өз «в өз «в «> өз өз өз өз «в «> өз е. «өз өз  ———— a —— ———————« 


Variables: 


plicnt - plot counter for hard copy displays 

Imagestr - string variable of the input data 

с0 - resolution level 0 data input 

cltocN  - corresponds to the c coeff level -] to -N 
dlltod1N - corresponds to the dl coeff's 

d21 to d2N - corresponds to the 12 соев 

d31to d3N - corresponds to the d3 coeff's 
Hh,Hv,Gh,Gv-  h and g coeff's in the horiz & vert direcuons 
HhHv,HhGv, 

GhHv,GhGv - 2-D decomposition masks 

lowest - lowest resoluuon level (a neg number) 
phsvct? — - phase vector for either the h or v direction 


эЧ эЧ ЭЧ эЧ эЧ ЭЧ ЭЧ АЧ АЧ эЧ aq ЭЧ ЭЧ GR OR АЧ АЧ АЧ эЧ GR OR GR GR GR a® af 


ооо ооо ооо ооо ооо ооо тосе 


% initialization of the counters and plots 


clg 

axis( normal’); 

hold off, 

рИсти= 0; 

chkputer 

clc 

ос EN. 
б.с есы ые ышаны ЕЕ 


* 2-D data input routine 


Imagestr =input(’ Enter the name of the image data: ','8"); 
с0 = eval(Imagestr); 
{il jl} =size(cO); 
if jl< =] | il< =] 
disp(’Bad Input Data, Restart the Program’) 
return 
end 
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€ 2-D wavelet coefficient input routine 
hcoefin 


% Branch for the 2-D multiphase case 


disp({’ — ']) 
q=input(’Do you desire the 2-D multiphase decomposition only? ’,’s’); 
ПЕ |4-- Y 

mp2d 


% Initialization for the decomposition routine 

% ofthe single selectable phase 

% 

% initializing Cout for the decomposition routine. Cout along 
% with Dlout, D2out and D3out are intermediate working 

% variables for the determining the coefficients for a 

% particular resolution level and phase. “phas??.m” are 

% one of for possible phase decomposition routines for the 
% two dimensional case. 


Cout=c0; 

subplot(211), mesh(cO0),utle('3-DPlot of the Data Array’) 
subplot(212),contour(cO),utle('ContourDisplay of the Data Array’) 
pause 

strplt 

clg 

clc 


dip(' — ']) 

q=[ Enter the number of resolution levels desired’ 
"Гог decomposition of the image: ae 

disp(q) 

lowest=input(’: > ^"); 

cle 

lowest =-abs(lowest); 


phsvcth = []; 
phsvctv 7 []; 


% The actual decomposition for the single 
% selectable phase case 


for lvl=-1:-l:lowest 


cle 
q— [' Enter the desired phase decomposition: ' 


'A) Horizontal Phase 0, Vertucal Phase 0' 
'B) Horizontal Phase 0, Vertical Phase 1” 
'C) Horizontal Phase 1, Verucal Phase 0' 
'D) Horizontal Phase 1, Vertical Phase 1 ' 
T 

disp" — Т) 

disp(q) 

phase =input(’: ',’s’); 


if phase==’A’| phase=='a’ 


phsvcth =(0 phsvcth]; 
phsvctv z [0 phsvctv]; 


Өз 


phas00 

elseif phase = =’B’| phase= =’b’ 
phsvcth =[0 phsvcth]; 
phsvctv =[1 phsvctv]; 
рһаз01 

elseif phase = 2'C'| phasez z'c' 
phsvcth zZ[1 phsvcth]; 
phsvctv z [0 phsvctv]; 
phas10 

elseif phase = =’°D’ | рћаве= = '' 
phsvcth=[1 phsvcth); 
phsvctv =[1 phsvetv); 
рһав11 

сіве 
disp(' Error") 
return 

end 


end 


% clean up for better memory utilization 
% This is primarily for 386 MATLAB 
clear wrkspc Dlout D2out D3 out Cout wrk a b q qq qqq ql 


% Energy check for the 2-D single selectable phase case 
enrg 2d 

% clean up for better memory utilization 

% This is primarily for 386 MATLAB 
clear Ek lvlabcj 


% Display of the phase sequence in reverse bit order as per the notes 
clc 

disp({’ "] 2) 

disp( Order of the Selected Phases from the highest to the lowest level’) 
disp(' 2) 

disp( Horizontal: Reading left to right corresponds to going from the’) 
disp(’ higher resolution to the lower resolution level”) 

disp({’ 7) 

disp(fliplr(phsvcth)) 

disp({’ Т) 

disp(' Verücal: Reading left to right corresponds to going from the") 
disp(’ higherer resolution to lower resolution level’) 

disp(U °)’) 

disp(fliplr(phsvctv)) 

dip(" Э 

disp(' < RET > to continue’) 

pause 


% Conduct the recomposition routine 
recon2d 


cle 
disp(' СТ) 
disp(' ROUTINE COMPLETE’) 


% The End of WAV2D.M 
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ыы со сооососеввеововьоовотеовоооосоосооооовњеоњо 


% HCOEFIN.m  2-D wavelet coefficient input 


% 

% Ву: LT J. E. Legaspi 

% Professor lam 

% US Naval Postgraduate School, Monterey CA 
% e-mail: legaspi@ece.nps.navy.mil 

% lam@ece.nps.navy.mil 

% Phone: (408) 646-3044/2772 


КРЕ 641916 1619191%1614 161% 9191010 91% 919196 19/4 6.4109 96 169 %/% 91910 0 9/0 % 9/06 10/6 6 % 0 9/96 0 96% 


% This routine obtains the desired wavelet h coefficients for both 
% the horizontal and vertical directions. The wavelets are assumed 


% separable. 

ТЕЛЯ Сао 

% Variables: 

% 

%  Hh,Hv,Gh,Gv- hand g coeff's in the horiz & vert directions 
€*  HhHv,HhGv, 

*  GhHv,GhGv - 2-D decomposition masks 

% wavelet? - string variable in the h or v directions 


% Obtaining the horizontal h coefficients 
disp([’ T) 
ql =[ Enter the number for the desired choice ” 
*for the wavelet in the HORIZONTAL direction:’ 
'(1) Haar h coefficients я 
'(2) Daubechies h coefficients : 
'(3) Own set of h coefficients ; 
I ' 
disp(q1) 
qqq 7 input(" "); 
ff — qqq--3. 
Hh= input(’ Enter your own set of H coefficients’); 
a=sum(Hh);b =Hh*Hh’; 
while a~ =] & b~ =.5 
disp(’ Your h coefficients are NOT compactly supported!’) 
Hh=input(’ Re-enter the h coefficients or ctrl-C to stop’); 
end 
waveleth =[’Own’); 
elseif qqq==2, 
99 = input('How many taps (2-10)? ^); 
Hh = daubdata(qq); 
waveleth=[’ Daubechies ’,num2str(qq),’-tap’}; 
else Hh=[.5 .5]; 
waveleth =[’Haar’}; 
аач=1; 


% Now for the vertical case 
disp([' 4 E) 
q 7 input('Do you desire the same wavelet in the vertical direction?’,’s’); 
if q==’Y’ | q == Evi 
Hv=Hh; 
wavelety = waveleth; 
else 
disp([’ D 
q! =[' Enter the number for the desired choice’ 
'for the wavelet in the VERTICAL direction: ’ 


'(1) Haar h coefficients 


'(2) Daubechies h coefficients : 

'(3) Own set of h coefficients i 

?, М: 
disp(q1) 
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qqq 7 input(" ); 
if qqq==3, 
Hv =input(’Enter your own set of H coefficients’); 
a =sum(Hh);b =Hh*Hh’; 
while a~ =1 & b~ =.5 
disp(’ Your h coefficients are NOT compactly supported!’) 
Hh =input(’ Re-enter the h coefficients or ctrl-C to stop’); 
end 
waveletv — ['Own']; 
elseif qqq==2, 
qq 7 input(' How many taps (2-10)? °); 
Hv =daubdata(qq); 
waveletv =[’ Daubechies ’,num2str(qq),’-tap’}; 
ele Hv=[.5.5]; 
waveletv — ['Haar']; 


адат 

end 
clc 
end 
Еа сас ee oe сс E 
Е.С E E ETE 
% Generating the g coefficients 
Gh =flipir(Hh); 


for i22:2:length(Hh) 
Gh(1,1) z-Gh(1,i); 

end 

Gv — fliplr(Hv); 

for i=2:2:length(Hv) 
Gv(1,1) =-Gv(1,1); 


end 

Qa c сесли сес сс ee 

REN ES Осы соске 

% Generating the 2-D masks for the decomposition 
Hv=Hv’; 

Gv =Gv’; 


% the 2° factor here reduces the number of multiplication? for 
* each lower level coefficients 

HhHv =2*Hv*Hh; 

HhGv =2*Gv*Hh; 

GhHv =2*Hv*Gh; 

GhGv =2*Gv*Gh; 

NHh =length(Hh); 

NHv =length(Hv); 

% End 


%%%%%009000060%09%9%40%440640060060000000еө909%%90%%%0%0%0%0%00000ФФОФФФОФФФФФе» 


* PHASO0.M 2-D phase 00 decomposition of the data array 


% 

% 

Х Ву: LT J. E. Legaspi 

% Professor Lam 

% US Naval Postgraduate School, Monterey CA 
% e-mail: legaspi@ece.nps.navy.mil 

% lam@ece.nps.navy.mil 

% Phone: (408) 646-3044/2772 


злой ө 


% This routine does a phase 00 decomposition given the input array - Cout. 
Ф horiz: Q vert: 0 


% 

* Variables: 

* Cout, Dlout-D3out - working variables for the coefficients 
* wrkspc - working array for the c coeff's 
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% Required m-files: 


% wav2d.m  - calling program 

% decomplt.m - plotting routine for the res level 
% 

ы еее не-е-е 
сіс 


[rows cols] 2 size(Cout); 
€ check if the row is odd, if it is, pad with one zero 
if rem(rows/2, floor{rows/2)) = = 
Cout=[Cout;zeros(1 ,cols)}; 
rows=rows+ 1; 
end 
%do the same with the columns 
if rem(cols/2,floor(cols/2)) = = 0 
Cout 7 [Cout zeros(rows, 1)]; 
cols 2 cols 4- 1; 
end 
% generate the workspace 
wrkspc 7 [zeros(NHv-1,cols -2* NHh-3);zeros(rows, NHh-1Cout zeros(rows, NHh-2);... 
zeros(NHv-1,cols - 2* NHh-3)]; 
clear Cout D1out D2out D3out 
% Now operate on the wrkspc with the decomposition masks 


a0; 
for row] =1:2:rows+NHv-2 
а=а+ 1; 
b=0; 
for coll =1:2:cola+NHh-2 
b=b+1; 


wrk = wrkspc(row1:row1 + NHv-1!,coll:coll + NHh-1); 
Cout(a,b) = sum(sum(HhHv .* wrk)); 
D1] out(a,b) =sum(sum(HhGv.* wrk)); 
D2out(a,b) =sum(sum(GhHv.* wrk)); 
D3 out{a,b) =sum(sum(GhGv.* wrk)); 
end 
end 
€ Store the working variables with the appropriate name 
eval(['c' ,num2su(abs(1v1)),' = Cout;’}) 
eval(['d1',num2stu(abs(1v1)),' 2Dl1out;']) 
eval(['d2'  num2str(abs(1v1)).' 2 D2out;']) 
eval(('d3' ,num2str(abs(1v1)),'  D3out;']) 
clg 
% Display the variables C, D1, D2 and D3 for the resolution level 
decomplt 


ооо ооо ео 


% РНА$10.М  2-D phase 01 decomposiUon of the data array 


E 


LT J. E. Legaspi 
Professor Lam 
US Naval Postgraduate School, Monterey CA 
e-mail: legaspiece.nps.navy.mil 
lam@ece.nps.navy.mil 
Phone: (408) 646-3044/2772 


ооо ооо ооо ооо 
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% This routine does a phase 10 decomposition given the input array - Cout. 
%  horiz: 1 vert: 0 


% 

% Variables: 

% Cout, Dlout-D3out - working variables for the coefficients 
% wrkspc - working array for the с сос Гв 


oF 


% Required m-files: 


% wav2d.m  - calling program 
* decomplt.m - plotting routine for the res level 
ЖК T E а M 


[rows cols] 2 size(Cout); 
% check if the row is even, if it is, pad with one zero 
if rem(rows/2,floor(rows/2)) = = 
Cout- [Cout;zeros(1 ,cols)]; 
rows=rows+ 1; 
end 
% odd column check 
if rem(cols/2,floor(cols/2)) ~ = 0 
Cout 7 [Cout zeros(rows, 1)]; 
cols 2 cols 1; 
end 
% generate the workspace 
wrkspc = [zeros( NHv-1 ,cols + 2* NHh-4) ;zeros(rows, NHh-2Cout zeros(rows, NHh-2),... 
zeros( NHv-2,cols+ 2*NHh-4)]; 
clear Cout Dlout D2out D3out 


а-0; 

for row] = 1:2:rows+ NHv-2 
а=а-+]; 
b=0; 
for coll =1:2:cols+NHh-2 
Ъ=Ъ+1; 


wrk 2 wrkspc(row1:row] -- NHv-1,coll:coll * NHh-1); 
Cout(a,b) 2 sum(sum(HhHv .* wrk)); 
D1out(a,b) 2 sum(sum(HhGv.* wrk)); 
D2out(a,b) - sum(sum(GhHv.* wrk)); 
D3out(a,b) 2sum(sum(GhGv.* wrk)); 
end 
end 

1/12 - abs(1v]); 

eval(['c',num2str(]1v12),' 2 Cout;']) 

eval(['d1',num2surt1vi2),'  D1out;']) 

eval(['d2' ,num2str(lvl2),' = D2out;’]) 

eval(('d3',num2sur(1vI2),' 2 D3out;']) 

clg 

decomplt 


© »»99999999999999999999999999999999999999999999999999999999999999999$5$ 


% PHASOI.M  2-D phase 01 decomposition of the data array 


* 

% 

% Ву: LT J. E. Legaspi 

% Professor Lam 

% US Naval Postgraduate School, Monterey CA 
% e-mail: legaspi@ece.nps.navy.mil 

% lam@ece.nps.navy.mil 

% Phone: (408) 646-3044/2772 


Mad ЕЕ ESI BO TB TE e I Ld LL EL LLL LL EI I I 


% This routine does a phase 01 decomposition given the input array - Cout. 
% horiz: 0 уеп: 1 


% 

% Variables: 

% Cout, Dlout-D3out - working variables for the coefficients 
% wrkspc - working array for the c сос (Гв 

OPT ee арасан л т 

% Required m-files: 

% wav2d.m  - calling program 

% decomplt.m- plotting routine for thc res level 
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сіс 
[rows cols] 2size(Cout); 
€ check if the row is odd, if it is, pad with one zero 
if rem(rows/2,floor{rows/2)) ~ = 
Cout 7 [Cout;zeros(1,cols)]; 
rows=rows+ 1; 
end 
% check if the columns are even 
if rem(cols/2,floor(cols/2)) = = 
Cout=[Cout zeros(rows,1)}; 
cols=cols +1; 
end 
% generate the workspace 
wrkspc =[zeros( NHv-2,cols +2* NHh-3);zeros(rows, NHh-1Cout zeros(rows, NHh-2);... 
zeros( NHv-2,cols + 2*NHh-3)}; 
clear Cout Dl out D2out D3out 
а=0; 
for row] =1:2:rows+NHv-2 
а=а+1; 
ь=0; 
for coll =1:2:cols+NHh-2 
b=b+1; 
wrk = wrkspc(row]:rowl + NHv-1,coll:coll +NHh-1); 
Cout(a,b) 2 sum(sum(HhHv .* wrk)); 
Dl1out(a,b) Zsum(sum(HhGv.* wrk)); 
D2out(a,b) =sum(sum(GhHv.* wrk)); 
D3out(a,b) 2 sum(sum(GhGv.* wrk)); 
end 
end 
1У12 — abs(lvl); 
eval(['c' , num2str(lv12),' 2 Cout; ']) 
eval(('d1',num2sut(lvD2),' 2 D1out;']) 
eval(['d2',num2str(1v12),' 2 D2out;']) 
eval(['d3' , num2su(lvD),' 2 D3out;']) 
clg 
decomplt 


адд 


% PHAS11.M 2-р рћһаѕе 11 decomposition of the data array 


% 

% 

% By: LT J. E. Legaspi 

% Professor Lam 

% US Naval Postgraduate School, Monterey CA 
% e-mail: legaspi@ece.nps.navy.mil 

% lam@ece.nps.navy.mil 

% Phone: (408) 646-3044/2772 


ооо ооо ссор ооо ооо хз 


% This routine does a phase 11 decomposition given the input array - Cout. 
%  honz:l vert: 1 


% 

% Variables: 

% Cout, Dlout-D3out - working variables for the coefficients 
% wrkspc - working array for the c coeff's 
%---.----------------.-.................................................... 
% Required m-files: 

% wav2d.m  - calling program 

% decomplt.m - plotting routinc for the res level 
% 

с 
сіс 


29 


[rows cols] 2 size(Cout); 


"e check of the row is odd, if it is, pad with one zero 
if rem(rows/2,floor(rows/2)) — = 0 
Cout =[Cout;zeros(1,cols)]; 
rows=rows+1; 
end 
%do the same with the columns 
if rem(cols/2, floor(cols/2)) ~ = 0 
Cout =[Cout zeros(rows, 1)]; 
cols =cols+ 1; 
end 
% generate the workspace 
wrkspc =[zeros(NHv-2,cols + 2* NHh-4);zeros(rows, NHh-2Cout zeros(rows, NHh-2);... 
zeros(NHv-2,cols -2* NHh-4)]; 
clear Cout D1out D2out D3out 


а=0; 
for row121:2:rows-- NHv-2 
а=а+ 1; 
Ь=0, 
for coll =1:2:cols+NHh-2 
b=b+1; 


wrk =wrkspc(row1:row1 + NHv-1,coll:coll +NHh-1); 
Cout(a,b) = sum(sum(HhHv .* wrk)); 
Dlout(a,b) =sum(sum(HhGv.* wrk)); 
D2out{a,b) =sum(sum(GhHv.* wrk)); 
D3 out(a,b) =sum(sum(GhGv.* wrk)); 
end 
end 

ІУІ? таБә(1МІ); 

eval(['c',num2strlvI2),' 2 Cout;']) 

eval(['d1',num2str(IvI2),' z D10out;']) 

eval(['d2' , num2str(Ivi2),' 2 D2out;']) 

eval(['d3' ,num2str(Ivl2),' 2 D3out;']) 

clg 

decomplt 


%%%9ч%Ф%0000000000009е0%00090900590е040940900000000%0е00000000е0000е40ФеФ000Фее0Ф%еФ0е 


% DECOMPLT.M displays the coefficients for a particular level 


% 

% Ву: LT J. E. Legaspi 

% Professor Lam 

% US Naval Postgraduate School, Monterey CA 
% e-mail: legaspi@ece.nps.navy.mil 

% lam@ece.nps.navy.mil 

% Phone: (408) 646-3044/2772 


( %86е%ғ80%в66566600%0 90 98 6 9 0 9 0 6960 8 9 9 9 0.0.0 6/6 6 9 е 0 6 9 69 9 0 9 6 4 9 9 0 0 0 6/0 0 9 6 9 9 9 6.6 6.6 9/6 е 


% Variables: as defined in WAV2D.M and PHAS??.M 


% 

% Required m-files: 

% WAV2D.M parent routine 

* PHAS??.M one of four possible calling routines 
% STRPLT.M metafile storage of selected plots 

y MM 

clg 


subplot(221),contour(Cout); xlabel('Cm"); 

utle([' wavelet horiz: ',wavelcth]); 

subplot(222),contour(D10ut); xlabel("D1m") 

utle([' wavelet vert: ',waveletv]) 

subplot(223),contour(D2out);title('D2m') 

xlabel(['Horiz phase: ',seteu(fliplr(phsvcth * 48)),' Vert phase: ',... 
sctsut(fliplr(phsvctv -48))]) 

subplot224),contour(D3out);utle('D3m) 
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xlabel([' Resoluuon Level ',num2sutlvl)]) 
pause 
вири 
сір 
subplot(221),mesh(Cout); xlabel(’Cm’) 
utle(['wavelet horiz: ',waveleth]); 
subplot(222),mesh(D1 out); xlabel(’D1m’) 
utle(['wavelet vert: ',waveletv]) 
subplot(223), mesh(D2out;Uütle('D2m") 
xlabel(['Honz phase: ',setsu(fliplr(phsvcth - 48)),' Vert phase: ',... 
setstr(fliplr(phsvctv +48))]) 
subplot(224), mesh(D3out);utle('D3m") 
xlabel([' Resolution Level ',num2sutlvl)]) 
pause 


suplt 


* END 


ИЛЫ Л УРУР зоо зе зәл» ези еее зул» 


% enrg2d.m Determining the Energy in each resolution level 

* Two Dimensional Case 

% 

% By: LT J E. Legaspi 

* Professor Lam 

% US Naval Postgraduate School, Monterey CA 

% e-mail: legaspi@ece.nps.navy.mil 

% lam(gece.nps.navy.muil 

% Рһопе: (408) 646-3044/2772 
ооо ооо ооо осо %% 
% For the single phase case, the energy is normalized to resolution 
% level 0. All the energies of the "c's" and the “d's” of a 

% particular level should add up to the energy of the "c's* of the 

% next higher level. The calling routine is “wav2d.m* and “strplt.m” 
% is a necessary m-file. 


РЕ ---................................... 
% Variables (other than defined in mp2d.m): 

% 

* E - energy array for [c,d1,d2,d3] for eaach level 

% esum - total energy of each level: [c(m) +d! +d2+d3]=c(m+1) 
ЕС --....................................... 


% Initialization of the energy sum for resolution level 0 
E=zeros(-lowest+ 1,4); 
E(1,1) zsum(sum(cO .^2)); 
ІМІ--1: 
k=l; 
% Calculate the energies for each resoluuon level 
while lvl > =lowest 
E(k - 1,1) Zsum(sum((eval(['c' , num2sut(abs(1vD))))). ^2)); 
E(k * 1,2) 2sum(sum((eval(['d]' , num2sut(abs(lv1))))) ^2)); 
E(k +1,3) 2 sum(sum((eval(['d2' , num2sut(abs(1vD)))) ^2)); 
E(k 1,4) 2 sum(sum((eval(['d3',num2sut(abs(lvD)])) ^2)); 
М1 =1У1-1; 
k=k+1, 
end 
E=E/E{1,1); % normalize the energy to level 0 
E=flipud(E); % top row is the lowest resolution 
* Display of the energies by coefficients 
clg 
k=[lowest-1 1 0 1.1]; 
axis(k);subplot(221),bar(lowest: 1:0, E(:,1)) 
xlabel{’c’) 
utle(['wavelet horiz: ',waveleth]) 
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axis(k);subplot(222), bar(lowest: 1:0, E(:,2)) 
xlabel('dl") 
utüe([' wavelet vert: ',waveletv]) 

axis(k);subplot(223) ,bar(lowest: ] :0, E(:,3)) 
utle('d2") 
xlabel('single phase case") 

axis(k);subplot(224) ,bar(lowest: 1:0, E(:,4)) 
tide('d3") 

pause 

вири 

сір 

сіс 

Е-Е” 

esum — sum(E); 

dip(' ']) 


аа = [' Verification of the energy in each resolution level ’ 
"Level Energy in Energy in C Difference’ 
' m level m-1 level m , 
j с+91+92+93 E 
и ']; 

disp(dd) 

] 7 -lowest; 


for i=1:-lowest 
a =esum(i);b =E(1,1+ 1);c =-j+1;d =abs(a-b); 
fprintf %2.0f %17.15Ғ %17.15f  ',c,a,b) 
fprintf(’ %17.15f\n’,d) 
j=j-1; 

end 

pause 

clc 


% End 


ALA EB HR AE IER HR EH ERIE IE IEEE IIE HEAL RE ES EAE AERA A IR ACE IEE 


% recon2d.m 2-D reconstruction algorithm for the single 

% selectable phase case 

% 

% 

% By: LT J. E. Legaspi 

Professor Lam 

US Naval Postgraduate School, Monterey CA 

e-mail: legaspi@ece.nps.navy.mil 
lam@ece.nps.navy.mil 

Phone: (408) 646-3044/2772 


ооо ооо ооо ооо = + т = = = = ® = = == т = = == о а а ә а 


Variables (other than defined in wav2d.m): 


HvHh,HvGh - reconstrution masks 

GvHh,GvGh 

whichl - phase decompositon for the particular level 
cwork - reconstructed c coefficient array for a level 


Required routines: 
wav2d.m 
strplt.m 
mmlt??.m -- ?? = 00,01,10,11 for the 
mask multiplication 


3% 33 эЧ эЧ 3% 33 33 33 33 35 38 34 33 мч эч 3S OS OX OR o 


РТТ ОТТО ТТТ ООВ 


б 
de 
о 


disp([' P 
q=input(’Do you desire to do the reconstruction (Y/N)? ’,’s’); 
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if q=="y’ СЕ =“ 
% The following four lines generate the reconstruction masks 
HvHh- fliplr(flipud(HhHv)); 
HvGh =flipIr(flipud(GhHv)); 
GvHh = fliplr(flipud(HhGv)); 
GvGh=flipIr(flipud(GhG v)); 
dip(' Т) 
q=input(’Do you desire to sce the masks? ','в'); 
ШЕ- x | q=='y’ 
% We will border the arrays with zeros for a better mesh plot 
[aa bb] 2size(HhH v); 
disp( ") 
disp HhHv decomposition") 
disp(HhHv) 
subplot(121) 
axis('square") 
mesh([zeros(1,bb 4 2);zeros(aa, 1)HhH v zeros(aa, 1);zeros(1,bb * 2)]) 
title( HhHv decomposition mask’) 
Фізр(” ”) 
disp('HvHh reconstruction!) 
disp(HvHh) 
subplot(122) 
axis('square') 
mesh([zeros(1,bb 4 2);zeros(aa, 1)H vHh zeros(aa, 1);zeros(1,bb - 2)]) 
utle(HvHh reconstruction mask") 
pause 
strplt 
сір 
disp(’ ’) 
disp(’HhGv decomposition’) 
disp(HhGv) 
subplot(121) 
axis("square’) 
mesh([zeros(1,bb + 2);zeros(aa, 1)HhGv zeros(aa,1);zeros(1,bb +2)]) 
title(’ HhGv decomposition mask’) 
disp(’ °) 
disp(’GvHh_ reconstruction’) 
disp(GvHh) 
subplot(122) 
axis('square") 
mesh([zeros(1 ,bb 4 2);zeros(aa, 1)GvHh zeros(aa,1);zeros(l, bb 2)]) 
tile('GvHh reconstruction mask") 
pause 
вирі: 
сір 
disp(" ") 
disp(’GhHv decomposition’) 
disp(GhHv) 
subplot(121) 
axis( square’) 
mesh([zeros(1 ,bb + 2);zeros(aa,1)GhHv zeros(aa,1);zeros(1,bb+2))) 
utle(’GhHv decomposition mask’) 
disp(’ ’) 
disp('HvGh reconstruction’) 
disp(H vGh) 
subplot(122) 
axis('square') 
mesh([zeros(1,bb + 2);zeros(aa, 1) HvGh zeros(aa, 1); zeros(1,bb 4 2)]) 
utle("HvGh reconstruction mask") 
pause 
strplt 
clg 
disp( ”) 
disp GhGv decomposition?) 
disp(GhGv) 
subplot(121) 
axis('square') 
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теяҺ((2егов(1,ЫЬ--2);гегов(аа,1)СОҺОу ?егов(аа,1);2егов(1,ЬЬ--2))) 
title(' GhGv decomposition mask") 
disp(" ") 
disp( 'GvGh reconstruction") 
disp(GvGh) 
subplot(122) 
axis('square") 
mesh([zeros(1,bb ^ 2);zeros(aa,1)GvGh zeros(aa,1);zeros(1,bb  2)]) 
utle(’GvGh reconstruction mask’) 
pause 
strplt 
clg 
end 
clear aa bb 
axis(’normal’) 
cle 
disp({’  ']’) 
dip NOW DOING THE RECONSTRUCTION’) 


% start at the lowest level 
% and calculate the c’s for the next level 
lowest — -abs(lowest); 
lvl=lowest; 
cwork —eval(('c' num2str(abs(lowest))]); 
аа =0; 
while lv1« 0 
аа =аа+ 1; 
(а Б|-віге(еуа((7с”,пшт2ви(аһө(1У1-- 1))))); 
d] zeval(('d1',num2str(abs(1v1))]); 
d2-— eval(['d2',num2str(abs(1v1))]); 
d3 — eval(('d3',num2str(abs(1v1))]); 
whichl 7 (num2str(phsvcth(a2)) ,num2stu(phsvctv(aa))]; 
eval(('tc 2 mmlit',which1,'(cwork, HvHhb);']) 
eval(('td] 2mmlt',whichl,'(d1,GvHh);']) 
eval(('td2 2 mmlt' ,whichl,'(d2, HvGh);']) 
eval(['td3 2 mmlt' ,whichl,'(d3,GvGh);']) 
clear cwork d1 d2 d3 
cwork ztc(1:8,1:b) *-td1(1:8,1:b) * td2(1:8,1:b) * td3(1:8,1:b); 
clear tc td] td2 td3 a b whichl 
clg 
axis('square") 
% plot the result of the c‘s for comparison 
subplot(121),contour(cwork); 
title(’ Reconstructed C array’) 
subplot(122),contour(eval([' c'  num2str(abs(lvl)- 1)])); 
title(('Actual C array']) 
xlabel(('Resolut;ion Level ',num2str(lv14- 1)]) 
pause 
strplt 
clg 
subplot(121),mesh(cwork); 
ütle(' Reconstructed C array') 
subplot(122) , mesh(eval(('c' , num2str(abs(lv1)- 1)])); 
ütle(('Actual C array']) 
xlabel([’ Resolution Level ',num2strlvl 4 1)]) 
pause 
strplt 
clg 
errr 2 cwork-eval((['c',num2str(abs(lvl 4 1))]); 
if max(max(abs(errr))) - =0 
mesh(abs(errr)) 
utle(('Absolute error from ',num2str(lvl),' to ',num2str(lvl* 1)]) 
xlabel(('Maximum Error Value: ',num2str(max(max(abs(errr))))]) 
pause 
strplt 
else 
clc 
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clg 
dispi' ОТ) 
dispC *** NO ERROR IN THE RECONSTRUCTION ***’) 
disp ТТ) 
disp NOW CALCULATING THE NEXT RECONSTRUCTION’) 
disp(U FROM LEVEL ‘num2str(ivi+1),’ TO ’,... 
num2sur(lvl 4 2)]) 

end 

ІМІ МІБ 1; 

end 

axis(’normal’) 

cle 

end 


function x =mml]t00(array ,mask) 


SSeS SSE CESSES SEES SS CESSES ESSERE SEES SEE SEE SESS S ESSE EEE ES 


*$ MMLTO0O.M | mask muluplicauon and summation for the 2-D phase 00 
% reconstruction 


% 

% Ву: LT J. E. Legaspi 

% Professor Lam 

% US Naval Postgraduate School, Monterey CA 

% e-mail: legaspi@ece.nps.navy.mil 

% lam@ece.nps.navy.mil 

% Phone: (408) 646-3044/2772 

Э © 20 1 2 2 eee SSS SSeS SSCS SESS SSS SESE SERS ESS SCS SSO S SCS SOS ES 
Ф ИНИИ Ты гоципе is only for PHASE-00 !!!!!!!!! 

% Required M - files: 

% wav2d.m -- grandparent routine 

% recon2d.m-- parent routine 

* ptzro2d.m-- puts n zeros rows & cols btwn coeff's 
Б ----------------.---------.------------------------.-- 


array —ptzro2d(array,1); 
% generate the workspace 
[NHv NHh] z size(mask); 
{row col] =size(array); 
array — [array zeros(row, NHh-2)];cols =col+ NHh-2; 
array — [array ; zeros(NH v-2,cols)]; 
а=0; 
for row1z1:row-1 
а=а+1; 
b=0; 
for coll 2 1:col-1 
b=b+1; 
wrk =array(row l: row] + NHv-1,coll:coll + NHh-1); 
x(a,b) =sum(sum(mask. *wrk)); 
end 
end 
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funcuon x - mmltlO(array, mask) 
omes такка а как кк кк кк мк к к кк т v 


* MMLTIO.M . mask multiplication and summation for the 2-D phase 10 
% reconstruction 


% 

% Ву: LT J. E. Legaspi 

% Professor Lam 

% US Naval Postgraduate School, Monterey CA 

% e-mail: legaspi@ece.nps.navy.mil 

% lam(gece.nps.navy.mil 

* Phone: (408) 646-3044/2772 
0o9999992992999999*39929999992*992999999229292999549€9929295999959*999$22292929959999929229295* 
% UH!!! This routine is only for PHASE-10 !!!!!!{!! 

% Required M=files: 

% wav2d.m -- grandparent routine 

% recon2d.m -- parent routine 

% ptzro2d.m-- puts n zeros rows & cols btwn coeff’s 
а 


array 7 ptzro2d(array,1); 
% generate the workspace 
[NHv NHh] 7 size(mask); 
[row col] — size(array); 
array — [zeros(row,1) array zeros(row, NHh-2); zeros(:NHv-2,col-- NHh-1)]; 
а=0; 
for row] =1:row-1 
а=а+1; 
b=0; 
for coll =1:col 
b=b+1; 
wrk = array(rowl:row1 + NHv-1,coll :coll + NHh-1); 
x(a,b) =sum(sum(mask. * wrk)); 
end 
end 


function x 2 mmlt01 (array , mask) 
qu ***9*99*9*5599*999299299*999992*55559929929929999992999999599395922395959599992929295999*9 


* MMLTO1.M mask multiplication and summation for the 2-D phase 01 
% reconstruction 


* 

* By: LT J. E. Legaspi 

% Professor Lam 

% US Naval Postgraduate School, Monterey CA 
* e-mail: legaspi@ece.nps.navy.mil 

* lam@ece. nps.navy.mil 

% Рһопе: (408) 646-3044/2772 


Сезе уен убу уубу өөө өөө уй 


Ф !'tt!t1!! This routine is only for PHASE-01 !!!!!t!!! 

% Required M = files: 

% wav2d.m -- grandparent routine 

% recon2d.m-- parent routine 

% ptzro2d.m-- puts n zeros rows & cols btwn coeff's 
еее. 

array = ptzro2d(array,1); 


% generate the workspace 
{NHv NHh] =size(mask); 
[row col] 2 size(array); 
array ={zeros(1,col+ NHh-2);arrayzeros(row, NHh-2);zeros(NHv-2,col+ NHh-2)]; 
а=0; 
for row] =1]:row 
а=а+]; 
b=0; 
for coll=1:col-1 
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b=b+1; 
wrk =array(row1:rowl + NHv-1,coll:coll + NHh-1); 
x(a,b) =sum(sum(mask. *wrk)); 
end 
end 


function x = mmlt1 1(array,mask) 
Re 7 7 CES SSST SS SS SSE SESS SSS S SESS SSS SS SESS SSE SSS SSH SSSSSS SSS SSSSS 


% MMLT11.M_ mask muluplication and summation for the 2-D phase 11 
% reconstruction 


% 

* By: LT J. E. Legaspi 

% Professor Lam 

% US Naval Postgraduate School, Monterey CA 

% e-mail: legaspi@ece.nps.navy.mil 

% lam@ece.nps.navy.mil 

% Phone: (408) 646-3044/2772 

АЛ ЛАлҚАҚААААЛАЛАЛЛАЛАЛААЛАЛЛАЛАЛАЛАЛАЛАЛАЛАААЛААЛАЛАЛЫ А. Л.А.АЛ.4.2,4.2.,4.2.2.2.4..2., 
% NN!!! This routine is only for PHASE-11 !!!!!!!!! 

% Required M - files: 

* wav2d.m -- grandparent routine 

% recon2d.m-- parent routine 

% ptzro2d.m-- puts n zeros rows & cols btwn coeff's 
РЕТТЕ ы ее... ea nacaas 


array =ptzro2d(array,1); 
% generate the workspace 
[NHv NHb] = size(mask); 
{row col] =size(array); 
array =[zeroa(1,col);array];rows =row + 1; 
array =[zeros(rows,1) array};cols=col+1; 
array — [array zeros(rows, NHh-2)];cols 2cols  NHh-2; 
array — [array ; zeros(NH v-2,cols)]; 
а-0; 
for row1-l:row 
a=a+l; 
b=0; 
for coll =1:col 
b=b+1; 
wrk =array(row1:row] + NHv-1,coll:coll +NHh-1); 
x(a,b) =sum(sum(mask.* wrk)); 
end 
end 
function x = ptzro2d(input,n) 


EE ЫК ЕКЕУ ЖУУИ ЕУ ЫК Көзөө өөө өөө ӨӨӨ 


%PTZRO2D.M puts "n" rows and columns of zeros after each 


% coefficient in the 2- D input array 

% 

* By: LT J. E. Legaspi 

* Professor Lam 

% US Naval Postgraduate School, Monterey CA 

* e-mail: legaspi@ece.nps.navy.mil 

% lam@ece.nps.navy.mil 

% Phone: (408) 646-3044/2772 
ооо соо охото се сое 
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% Required routines: 


% wav2d.m _ -- great-grandparent routine 
% recon2d.m_ -- grand parent routine 
% mmlt??.m — -- parent rouüne, ?? 200,01,10,11 
“лс с > er 
(г с|-віге(іпрш); 
x0 z[]; 
for j21:c 
х0 = [х0 input(: ,j) zeros(r,n)); 
end 
[r c] 2 size(x0); 
х= []; 
for j2l:r 
x z[x;xO(j, :);zeros(n,c)]; 
end 


mL SiS See Sse I DI c ВЕ P p Rp кан eee eee ee 


% MP2D.M 2-D Multü-phase decomposition 


% 

% Ву: LT J. E. Legaspi 

% Professor Lam 

% US Naval Postgraduate School, Monterey CA 
% e-mail: legaspi@ece.nps.navy.mil 

% lam@ece.nps.navy.mil 

% Phone: (408) 646-3044/2772 


ооо оао 


* MP2D conducts a multiple phase decomposition routine for the 2-D 
% case. The following m-files are required: 


% 

$ отау24лп - саШпр гошіле 

% | mpdcmp2d.m - multiple phase decomp for 1 resolution level 
%  strplt.m  - metafile storage of selected plots 

*  enrg2dmp.m - muluple phase energy check 
Пе 
% Variables are the same as the calling program 

% 

о тағ 
о ras) ce с 

% Initialization of the routine 

cle 

сір 


subplot(211), mesh(cO),Utle('3-DPlot of the Image Array’) 
subplot(212),contour(cO),utle(’ContourDisplay of the Image") 
pause 

strplt 

clg 

clg 

axis( normal") 

hold off 

dispi' — ']) 

lowest — input(' Enter the lowest resolution level: '); 
lowest — -abs(lowest); 

m —0, 


% Calculating the Coefficients for each resolution level 

while m> =lowest+ 1 
eval({’c’ ,num2str(-m+ 1),’=mpdcmp2d(c’ ,num2str(abs(m)),’, HhHv,m);’)) 
eval([’d1',num2str(-m + 1),’=mpdcmp2d(c’ ,num2str(abs(m)),’, HhGv,m);’)) 
eval(('d2' ,num2str(-m + 1),’=mpdcmp2d(c’ ,num2str(abs(m)),’, GhHv,m);']) 
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eval(['d3' ,num2str(-m * 1),' - mpdcmp2d(c' , num2su(abs(m)).' ,GhGv,m);']) 
% Display of the coefficients 
subplot(221), eval(['mesh(c',num2str(abs(m) * 1),')']) 
xlabel([’c’]) 
title([’wavelet horiz: ’ ,waveleth]) 
subplot(222),eval([’mesh(d1’ ,num2str(abs(m) + 1),")’]) 
xlabel(['d1']) 
utle(['wavelet vert: ',waveletv]) 
subplot(223),eval(['mesh(d2' ,num2str(abs(m) * 1).)']) 
utle(['d2']) 
xlabel('muluphase case") 
subplot(224), eval([('mesh(d3' ,num2str(abs(m) 4 1),')'] 
utle(['d3']) 
xlabel([' Resolution level ',num2str(m-1)]) 
pause 
strplt 
clg 
subplot(221), eval(['contour(c' , num2sur(abs(m) * 1),)']) 
xlabel(['c']) 
utle([' wavelet horiz: ',waveleth]) 
subplot(222), eval(['contour(d1' , num2str(abs(m) * 1), ')'] 
xlabel(['d1']) 
utle([' wavelet vert: ',waveletv]) 
subplot(223), eval(['contour(d2' ,num2str(abs(m) t 1),)'] 
utle(['d2']) 
xlabel(’multiphase case’) 
subplot(224),eval([’contour(d3’ ,num2str(abs(m) + 1).’)’]) 
ше(Г437) 
xlabel(['Resoluuon level ',num2str(m-1)]) 
pause 
strplt 
clg 


function x =mpdcmp2d(array ,mask,m) 

EE Ue sesse*tesS92999$999999909099999099990999299090999299909999999990999999 
% MPDCMP2D.M 2-D Multi-phase decomposition decomposition function 

% m-file 

% 

% Ву: LT J. E. Legaspi 

% Professor Lam 

% US Naval Postgraduate School, Monterey CA 

% e-mail: legaspi@ece.nps.navy.mil 

% lam@ece.nps.navy.mil 

% Phone: (408) 646-3044/2772 

EN LUIS T "aa aaonee"os0o99999990999999999999999999099999090999999 
% MPDCMP2D is the actual function m-file that zero pads the input 

unage array given the size of the mask operator and the 

resolution level m. "m" is the current resolution level 

and is a negative number. 


mp2d.m  - calling rouune 


a& а з al 94 384 a 


еч ен ән нен қ ңе өн че че чече че че че — а ча ар ча а ча н ча Ч на на на ча Ч ча о ца а р ча а в В в Че а momo TE eT eet 
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[NHv NHh] -size(mask); 
[a b] =size(array); 
x =zeros(a+2*(-m)*(NHv-]),b+2°(-m)*(NHb-1)); 
for k 1: NHh 
for 12 1: NHv 
%add the rows of zeros 
work 7 [zeros(2^(-m)* (1-1), b);array;zeros((NH v-1)*2" (-m).b)]; 
[a b] 2 size(work); 
%add the cols of zeros 
work] =[zeros(a,2“(-m)*(k-1))work zeros(a,(NHh-k)*2°(-m))]; 
work2 2mask(NHv-l- 1, NHh-k * 1)*workl ; 
x=x+work?2; 
€ clear work work1 work2 
end 


end 


CAE ITEM HM BE IR IOI D EC ака кт ктк та кк кек ал титан кок аа 


% enrg2dmp.m Determining the energy in cach resolution level 


% multiphase case 

% 

* By: LT J. E. Legaspi 

% Professor Lam 

% US Naval Postgraduate School, Monterey CA 
% e-mail: legaspi@ece.nps.navy.mil 

% lam(gece.nps.navy.mil 

% Phone: (408) 646-3044/2772 


Сезнеке е Ееее коео еы ороо ресе еее те етее RR I 


% In this MULTIPHASE case, the lower level will have 

% four possible phases, thus its energy will be 4 times as high as 
% the next higher level so we will divide the energy of the lower 
* level and the value should equal the energy of the c’s of the 

% next higher level. 


% Variables (other than defined in mp2d.m): 

% 

*$ E - energy array for [c,d1,d2,d3] for eaach level 

* esum - total energy of each level: [c(m)+d1 +d2+d3]/4 =c(m+1) 


E-zeros(-lowest 4 1,4); 

% Initialization of the energy for resolution level 0 

E(1 ,1) =sum(sum(c0 .*2)); 

lvl=-1; 

k-1; 

% Calculate the energies for each resolution level 

while lv] > =lowest 
E(k* 1,1) Z saum(sum((eval([' c' o num2str(abs(1v1))])). ^ 2*4* (1v1))); 
Е(К + 1,2) Zsum(sum((eval(['d1' , num2str(abs(1v1))])) ^2*4 ^(1v))); 
E(k * 1,3) Z sum(sum((eval(['d2' , num2str(abs(1v]))])) ^2*4^(14D)); 
E(k +1,4) Z eum(sum((eval(['d3' ,num2str(abs(1v1))])) ^2*4 ^ (1v1))); 
ІУІ-МІ-1; 
k=k+1; 

end 

E=E/E(1,1); % normalize the energy to level 0 

E=flipud(E); % top row is lowest resolution level 

% Display of the energies by coefficients 

clg 

k=[lowest-1 1 0 1.1]; 

axis(k);subplot(221), bar{lowest: 1:0, E(: ,1)) 

xlabel('c') 
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uue([' wavelet horiz: ',waveleth]) 
axis(k);subplot(222),bar(lowest:1:0,E(:,2)) 
хізБе (741?) 
utle([' wavelet vert: ',waveletv]) 
axis(k);subplot(223),bar(lowest:1:0, E(:,3)) 
utle('d2") 
xlabel('multüiphase case’) 
axis(k);subplot(224),bar(lowest: 1:0, E(:,4)) 
ütle('d3") 
pause 
strpit 
сір 
сіс 
Е-Е” 
esum —sum(E); 
esum —sum(E); 


disp °T) 


dd=[’ Verification of the energy in each resolution leve) — ' 
"Level Energy in Energy in C Difference’ 
"m level m-1 = level m 
li c+d1+d2+d3 қ 
ИА J 

disp(dd) 

j7 lowest; 


for i=1:-lowest 
a=esum(i);b= E(1 ,i+ 1);d =abs(a-b),;c=j+1; 
Ергіп!(” %2.00: %17.15f %17.15f ’,c,a,b) 
fpnntf(’ %17.15f\n’,d) 
]=)+1; 

епа 

pause 

clc 


% Ела 


ЕЕЕ 1 57%%%%59%8%%8%%%9%%%99%%9%0%9%0%99%%%%%9099%%90%0%400%0000%0000000000ө 


% IDATA.M Sample Image data for testing the routines 


% 

% By: LT J. E. Legaspi 

% Professor Lam 

% US Naval Postgraduate School, Monterey CA 
% e-mail: legaspi@ece.nps.navy.mil 

% lam@ece.nps.navy.mil 

% Phone: (408) 646-3044/2772 


E + oe SS SSIES CSHSESSSESS LSS SS SS SS VSS SSSSSS SSS SS SS SSS HSS SS SS SSS se Se 
K 2 menu('Select image','Box1','Box2','Box3','B/W','2B/W','8B/W',... 
’Х’, 'Нехароп’); 


Е = = 

* first image is a box 

im —2eros(100,100); 
im(20:21,20:80) 2 ones(2,61); 
im(79:80,20:80) 2 0nes(2,61); 
im(22:78,20:21) =ones(57,2); 
im(22:78 ,79:80) =ones(57,2); 


elseif K = = 

im =zeros(100,100); 
im(20,20:79) 2 ones(1,60); 
im(80,20:79) 2 ones(1,60); 
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im(20: 79,20) — ones(60,1); 
1m(20:79,80) = ones(60, 1); 


elseif K== 

im =zeros(100, 100); 
im(1:2,1:61) =ones(2,61); 
1m(60:61,1:61) 2 ones(2,61); 
im(3:59,1:2) z ones(57,2); 
im(3:59,60:61) =ones(57,2); 


elseif K= =4 

% second image is a 1 white and 1 black rectangle 
im =zeros(100,100); 

im(:,51:100) 20nes(100,50); 


elseif K= =5 

% third image are two black boxes and two white boxes 
un =zeros(100,100); 

im(1:50,1:50) =ones(50,50); 
im(51:100,51:100) = ones(50,50); 


elseif K= = 

% fourth image are 8 black and 8 white boxes 

im =zeros( 100,100); 

im(1:25,:) 2 [ones(25,25)zeros(25,25) ones(25,25) zeros(25,25)]; 
im(26:50,:) 2 [zeros(25,25) ones(25,25) zeros(25,25) ones(25,25)]; 
im(51:75,:) Z(ones(25,25)zeros(25,25) ones(25,25) zeros(25,25)]; 
im(76:100,:) =[2егоз(25,25)опез(25,25) 2егоз(25,25) ones(25,25)]; 


elseif K= = 

% fifth image are 2 45 degree diagonals 
im — eye(100,100); 

im —sign(im + flipud(im)); 


elseif K — — 
% sixth image is a hexagon 
for rowz 1:20 
в=2*го\/-1; 
im6(row,n:n * 2) 2 ones(1,3); 
end 
[a b] 2 size(im6); 
row2 —ceil(sqrt(a^2 4 b^2)); 
im6a —zeros(row2,b); 
im6a(:,b-1:b) Zones(row2,2); 
rim 7 [im6;im6a; fliplr(im6)]; 
lim = flipIr(rim); 
im6 — [lim rim]; 
clear row row2 im6a lim rim 
[a b] 7 size(im6); 
im6a =zeros(100,100); 
ва = Поог((100-а)/2); 
bb = Поог((100-Ь)/2); 
im6a(aa + 1:aa+a,bb+ 1:bb+b)=im6; 
im =im6a; 
clear im6 ипба а aa b bb K n 
end 


ІШ) 


C. ROUTINES COMMON TO BOTH DIMENSIONS 


Сезе зз зз өөө өөө ӨӨӨ Ө 


% chkputer.m checks for the proper operating system 


% 

% Ву: LT J. E. Legaspi 

* Professor Lam 

* US Naval Postgraduate School, Monterey CA 
% e-mail: legaspi@ece.nps.navy.mil 

% lam@ece.nps.navy.mil 

% Phone: (408) 646-3044/2772 


КЛ ЕУ Рекетке рекке» 


% 

% This routine reads the variable "cmptr" to see if it can run in 

* the current system. At this time, the algorithm will run on a 

% PC386, UNIX, or SUN system. If the operating system is a VAX- 
% based, the routine will not be as reliable since the VAX MATLAB 
% does not recognize IEEE arithmetic. The routine also deletes all 

% other leg*.met files prior to running. 


[cmptr mxsz] 2 computer; 
a =length(cmptr); 
Cic; 


if emptr(1)==’P’& cmptr(2) = 2'C'&cmpu(3) 2 =’3’, 
'del leg*. met; 
elseif cmptr(1) = 2'S'&cmptr(2) - -'U'&cmpu(3)2 2'N', 
Irm leg*.met; 
elseif cmpu(a-3) Z Z'U'&cmptu(a) 2 2'X', 
'rm leg*. met; 
else 
clc 
dip[' °T) 
qi z['*****Errorin input******' 
"Hardware cannot support ' 
“the algorithm... PLS A 
' Restart the Program =’); 
disp(q1) 
retum 
end 


О $ ооо ооо ооо ооо ооо тоне 


% вири. т 

% 

% Store Plot M-file for wavelet decomposition routine 
% 

* By: LT J. E. Legaspi 

% Professor Lam 

% US Naval Postgraduate School, Monterey CA 
% e-mail: legaspi@ece.nps.navy.mil 

% lam@ece.nps.navy.mil 

% Phone: (408) 646-3044/2772 


7.122224414444444444144424444444444444222444,,,...0.202002002022020.022222.2.220.02.02020202020:42.:0 


% 

% This program stores the plot as a metafile with the prefix "leg" 

% followed by a number starting from 0 on upwards and ending with 
% the extension ".met" You must be familiar with GPP to get hard- 
% copies of the plots. "pltcnt" is the only variable necessary for 

% the proper indexing of the leg*.met files. Ensure that the 
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% calling program does not redefine the variable! 


D EES o9 occ Аты 
disp({’ СҮ) 
q5 2 input('Store the Plot (Y/N)? ','8); 

if q5 = = '\’ | 45 =='y', 


eval([’meta leg’ num2str(plicnt)]) 
plicnt=pltcnt+ 1, 
end 


зз» өө лөө» өөө өөө өөө 


% daubdata.m 

% 

% By: CAPT/USMC Kalmbach, M. 

% Professor Lam 

% US Naval Postgraduate School, Monterey CA 
% e-mail: lam@ece.nps.navy.mil 

% Phone: (408) 646-3044 

% 

% 


funcuon[hcoeff] 2 daubdata(x) 
% This function retums the proper "b" coefficients for the specified order 
% Daubechies wavelet. The input value of "x" is the specified order. 
Í x == 
ВсоеЙ = [0.482962913145;0.836516303738;0.224143868042;-0.129409522551]; 
elseif x == 
hcoeff = [0.332670552950;0.806891 509311 ;0.459877502118;-0.135011020010; 
-0.085441273882;0.035226291882]; 
elseif x == 
hcoeff = (0.230377813309;0.714846570553;0.630880767930; -0.027983 769417; 
-0.187034811719;0.030841381836;0.032883011667;-0.010597401785); 
elseif x == 
hcoeff — [0.160102397974;,0.603829269797;,0.724308528438;0.138428145901; 
-0.242294887066;-0.032244869585;0.077571493840;-0.006241490213; 
-0.012580751999;0.003335725285]; 
elseif x == 
hcoeff — [0.111540743350;0.494623890398;,0.751133908021;0.315250351 709; 
-0.226264693965;-0.129766867567;0.097501605587;0.027522865530; 
-0.031582039318;0.000553842201;0.004777257511;-0.001077301085); 
elseif x == 
hcoeff — [0.077852054085;,0.396539319482;0.729132090846;0.469782287405; 
-0.143906003929;-0.224036184994:0.071309219267;0.080612609151; 
-0.038029936935;-0.016574541631;0.012550998556;0.000429577973; 
-0.001801640704;0.000353713800]; 
elseif x == 
hcoeff = [0.054415842243;0.312871590914;0.675630736297;0. 585354683654; 
-0.015829105256,-0.2840 15542962;0.000472484574;0. 128747426620; 
-0.017369301002;-0.044088253931;0.013981027917;0.008746094047; 
-0.004870352993;-0.000391740373;0.000675449406;-0.000117476784); 
elseif x == 
hcoeff — [0.038077947364;,0.243834674613;0.604823123690;0.657288078051; 
0.133197385825;-0.293273783279;-0.096840783223;0. 148540749338; 
0.030725681479;-0.067632829061;0.000250947115;0.022361662124; 
-0.004723204758;-0.004281503682;0.001847646883;0.000230385764; 
-0.000251963189;0.000039347320]; 
elseif x == 10 
hcoeff = [(0.026670057901;0. 188176800078 ;0.527201 188932;0.688459039454; 
0.281172343661 ,-0.249846424327;-0.195946274377;0. 127369340336; 
0.093057364604;-0.071394147166;-0.029457536822;0.033212674059; 
0.003606553567;-0.010733175483;0.001395351747;0.001992405295; 
-0.000685856695;-0.000116466855;0.000093588670;-0.000013264203}; 
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сізе 


hcoeff = (J; 

break 
end 
hcoeff = hcoeff'/sqrt(2); % normalize the “h” vector 
retum 


SSSSSHS SS SHSSSSSSSSSHS SS SHS SSS SSS SSS SHS SSS SSS SS HSS HSS SSH SSS SSH SHSSSS SKS HSE KSE SEES 


#! /bin/csh -f 

ф 

f**"**9*9*9*99*9*»559*5959559*** 5559 *5*5€5*59**9**99*555*9**€**55***t*sssstsss*sstttsssstut 
# metal4 C-Shell Routine 

й 


# metal4 is a command executable file for the generation of 
# the leg*.met metafiles into postscript graphic files for the 
# SPARC printer number 14 on the fourth floor of Spanagel 431. 
# This routine is system dependent and is guaranteed to work 
# on the UNIX based system of the ECE Department at the Naval 
# Postgraduate School. 
ооо азота оо 
# 
ls -1 leg*.met 2 tmmmpl 
echo "4! /bin/csh -f * > tmmmp2 
awk -F. '(printf "gpp * $1 ° -dps -ol \n” ;\ 
printf “Ipr14 * $1 °.ps \n";\ 
printf “rm ° $1 ".ps i")' tnmmpl » » tmmmp2 
chmod 0700 tnmmp2 
tmmmp2 
rm tmmmp1 tnmmp2 
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